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We study the nonperturbative dynamics of the standard model (SM) after inflation, in the regime 
where the SM is decoupled from (or weakly coupled to) the inflationary sector. We use classical 
lattice simulations in an expanding box in (3+1) dimensions, modeling the SM gauge interactions 
with both global and Abelian-Higgs analogue scenarios. We consider different postinflationary 
expansion rates. During inflation, the Higgs forms a condensate, which starts oscillating soon after 
inflation ends. Via nonperturbative effects, the oscillations lead to a fast decay of the Higgs into 
the SM species, transferring most of the energy into Z and bosons. All species are initially 
excited far away from equilibrium, but their interactions lead them into a stationary stage, with 
exact equipartition among the different energy components. From there on the system eventually 
reaches equilibrium. We have characterized in detail, in the different expansion histories considered, 
the evolution of the Higgs and of its dominant decay products, until equipartition is established. 

We provide a useful mapping between simulations with different parameters, from which we derive 
a master formula for the Higgs decay time as a function of the coupling constants, Higgs initial 
amplitude and postinflationary expansion rate. 

I. INTRODUCTION 

Inflation, an early period of accelerated expansion, is 
the leading framework to explain the initial conditions of 
the Universe. The concrete particle physics realization 
of inflation has eluded any clear identification so far, so 
the inflationary dynamics is often described in terms of a 
scalar field, the inflaton, with a vacuum-like energy den¬ 
sity. Furthermore, the confirmed discovery of the stan¬ 
dard model (SM) Higgs in the Large Hadron Collider 
(LHC) [HU] has initiated the quest for understanding the 
cosmological implications of the Higgs, and in particular, 
its possible role during and after inflation. Intriguingly, 
the SM Higgs could play the role of the inflaton, if a non- 
minimal coupling to gravity is introduced, appropriately 
fixed to fit the observed amplitude of the cosmological 
perturbations [3]. This model, known as Higgs-inflation, 
constitutes undoubtedly one of the most attractive and 
economical scenarios for realizing inflation, though it is 
not free of criticism mil]; see, however. Ref. [6]. 

In this paper, we will rather explore a different route 
for the possible role of the Higgs during and after infla¬ 
tion. We will merely assume that inflation was driven by 
a very slowly evolving energy density, without specifying 
the nature of the field responsible for it. Inflation can 
then be seen effectively as a quasi-de Sitter background 
with a slowly changing Hubble rate. We will assume that 
the SM Higgs is not coupled directly to the inflationary 
sector [THin]- Under these circumstances, the Higgs be¬ 
haves simply as a spectator field living in a (quasi-)de Sit¬ 
ter background, with the effective potential of the Higgs 
ultimately dictating its behavior. Let us note that even 
if there is no direct coupling, it is likely that effective 
operators will connect the Higgs with the inflaton, via 
some possible mediator field(s). Moreover, the need to 


reheat the Universe after inflation requires somehow the 
presence of such coupling, though there is no particu¬ 
lar constraint on this. As we will see, the Higgs decays 
very fast after inflation into all SM species, so one can 
safely assume that the effect of an inflaton-Higgs cou¬ 
pling is negligible, unless that coupling is significantly 
large. Therefore, even if such coupling is present, we will 
consider it weak enough so that any Higgs-inflaton inter¬ 
action does not affect the dynamics of the latter during 
or after inflation. 

The improved renormalized Higgs potential has been 
computed at next-to-next-to-leading order [TTl [T^ . It is 
characterized by the running of the Higgs self-coupling 
A(/r), which decreases with energy dX/dpi < 0, and be¬ 
comes negative above a certain critical scale ptQ, A(/i > 
To) < 0. Equivalently, the effective potential develops 
a barrier at large field amplitudes, reaching a maximum 
height at some scale t+ < To, so that at higher energies 
T > T+ the effective potential goes down, crosses zero 
aX T = To and becomes rapidly negative, possibly reach¬ 
ing a (negative) minimum at some scale T- ^ To- This 
can be seen in Fig. [l] These scales depend sensitively 
on the Higgs mass ttih, the strong coupling constant Og, 
and especially on the top Yukawa coupling yt. For the 
SM central values, Og = 0.1184, niH = 125.5 GeV, and 
the most recent measurement of the top quark mass by 
GMS, rut = 172.38 GeV [13], one finds /i_|_ ~ 2 x 10^^ 
GeV and /ig ~ 3 x 10^^ GeV. If one takes the world av¬ 
erage top quark mass rrit = 173.36 GeV El: then /r+, yo 
are reduced by a factor ~ 1/30. However, by considering 
a value of yt merely 2-3 sigma smaller than its central 
one, we can push the critical scales to t+,To ^ 5 x 10^® 
GeV. Besides, minimal additions to the SM such as a 
scalar singlet coupled to the Higgs [T5| SS] , or even a 
small nonminimal coupling of the Higgs to gravity 
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FIG. 1: Improved renormalized Higgs potential at next-to- 
next-to-leading order (red continuous line) computed for Qs = 
0.1184, niH = 125.5 GeV, and mt = 171.2 GeV (< 2a below 
CMS central value). Also shown, for comparison, the function 
(blue dashed line), where A+ = A(/i+) ~ 3 x 10“®. 


can also modify the running of A(/x) and stabilize the ef¬ 
fective potential. In such a case, the Higgs self-coupling 
may remain always positive A(/r) > 0. 

We will consider that the Higgs amplitude during in¬ 
flation remains always in the ’safe’ side of the effective 
potential, where A(/r) is positive. This can be guaranteed 
if is sufficiently large (compared to the inflationary 
scale), or alternatively, if beyond-the-SM physics stabi¬ 
lizes the potential at high energies. With these considera¬ 
tions, the Higgs fluctuates during inflation, like any light 
degree of freedom. The fluctuations then pile up at super- 
Hubble scales, creating a condensate [mun]. The am¬ 
plitude of the Higgs condensate, however, will not grow 
unbounded with the numbers of e-folds, as it happens in 
the case of a massless free field. On the contrary, the 
Higgs self-interactions provide an effective (sub-Hubble) 
mass to the fluctuations, which eventually saturates the 
growth of the condensate amplitude m- The distribu¬ 
tion of the Higgs amplitude at super-Hubble scales enters 
very fast, within a few e-folds, into a self-similar regime, 
which continues until the end of inflation. The Higgs 
condensate acquires this way a fixed physical correlation 
length (exponentially larger than the Hubble radius) and 
a fairly large amplitude. This will set up the initial con¬ 
dition for the behavior of the Higgs after inflation. 

Notice, however, that our analysis will not really de¬ 
pend on the condition A(/x) > 0 during inflation. The 
possible implication of the Higgs self-coupling becoming 
negative, A(/r) < 0, during inflation has indeed been an¬ 
alyzed in detail in Refs. [71EIH28]. In this case, if the 
scale of inflation is sufficiently high, the second minimum 
can be reached and anti-de Sitter bubbles are formed 
during inflation. The consequence of this does not need 
to be catastrophic, but rather indicative that either the 
condition that our Universe is in the electroweak (EW) 
vacuum is something very special (very improbable), or 


either some new physics beyond the SM is necessary to 
stabilize the Higgs potential. The crucial ingredient for 
our analysis is, therefore, not that the Higgs self-coupling 
A(/r) remains positive during inflation, but the fact that 
the Higgs develops a vacuum expectation value (VEV) 
during inflation much larger than the electroweak (EW) 
scale ~ 0(10^) GeV. The way such condensate is at¬ 
tained is mostly irrelevant. The case A(/r) > 0 all through 
inflation provides simply a reference case, where the for¬ 
mation of a large VEV during inflation is unavoidable, 
and its typical amplitude can be easily calculated. 

In this paper we investigate in detail the Higgs’s de¬ 
cay into its most energetically dominant decay products, 
the SM electroweak gauge bosons, during the immediate 
stages following the end of inflation. Our work represents 
a complementary analysis to that of Enqvist et al. mm, 
where analytical techniques were employed to study the 
same problem. We use instead lattice simulations in an 
expanding box in (3-1-1) dimensions, modeling the SM in¬ 
teractions with global and Abelian-Higgs setups, which 
go beyond the assumptions behind any analytical cal¬ 
culation. Besides this, we also consider different Higgs 
initial amplitudes and postinflationary expansion rates. 

The paper is organized in such a way that we increase 
progressively the complexity of the different approaches 
used to describe the dynamics of the system, approxi¬ 
mating the structure of the SM interactions better and 
better at each new step. In Section]^ we first present a 
brief analysis of the behavior of the Higgs after inflation, 
ignoring its coupling to the rest of the SM species. In Sec¬ 
tion m we switch on the coupling to the SM helds, but 
ignore the gauge nature of the interactions. We obtain 
analytical estimates for a later comparison with numer¬ 
ical simulations. In Section |IV| we present the first set 
of lattice simulations, where we follow the Higgs and its 
decay products, yet under the assumption that the gauge 
nature of the SM interactions can be neglected. In Sec¬ 
tion |Vj we finally incorporate gauge interactions into the 
simulations, by modeling the SM with an Abelian-Higgs 
setup. Although this is just an approximation to the 
gauge structure of the SM, the non-Abelian nonlineari¬ 
ties can arguably be neglected. Therefore, the outcome 
of these simulations represents the most precise calcu¬ 
lation of the dynamics of the SM after inflation, fully 
incorporating the nonlinear and nonperturbative effects 
of the SM, while considering the gauge nature of its inter¬ 
actions. In Section [VT| we present a useful mapping be¬ 
tween simulations with different parameters, from which 
we obtain a characterization of the Higgs decay width as 
a function of the coupling constants, initial Higgs ampli¬ 
tude, and postinflationary expansion rate. In Section [VH| 
we conclude and discuss some of the possible cosmologi¬ 
cal implications of our results. 

All through the text H = c = 1, and nip ~ 2.44 x 
10^® GeV is the reduced Planck mass. We take the flat 
Friedmann-Robertson-Walker (FRW) line element ds^ = 
a?(t){—dt^+dx'^dx'^) for the background metric, with a{t) 
the scale factor and t the conformal time. 
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II. HIGGS OSGILLATIONS AFTER INFLATION 


Let us characterize inflation as a de Sitter period with 
Hubble rate iJ* MeW) where Mew 0(10^) GeV 
is the EW scale. In reality, we know that inflation can¬ 
not be exactly a de Sitter background, since inflation 
must end after a finite number of e-folds. The curva¬ 
ture perturbation spectral index Ug = 0.968 ± 0.006 
constrained by Planck to be smaller than unity at more 
than 7 sigma, is actually interpreted as an indication of 
the quasi-de Sitter nature of Inflation. For our purposes, 
however, the distinction between de Sitter and quasi-de 
Sitter is irrelevant. 

With a gauge transformation, the SM Higgs doublet 
can be parametrized in the unitary gauge by a single 
scalar real degree of freedom, $ = <^1^/2. The renormal¬ 
ized improved potential for large-held amplitudes Iv^l ^ 
Mew is just given by the quartic part 

V(v) = (1) 

with A(/r) the Higgs self-coupling at the renormalization 
scale fj, = ip. Radiative corrections to the potential are 
encoded in the running of A(/r), which to date has been 
computed to three loops when the Higgs is minimally 
coupled to gravity mmi]. 

We ignore the nature of the sector responsible for inha- 
tion, so a priori there is no need for the Higgs to be cou¬ 
pled directly^ to the inhationary sector. We will just con¬ 
sider that the Higgs held simply ’lives’ on the de Sitter 
background, playing no dynamical role during inhation, 
and behaving simply as a spectator held [7HS]. As men¬ 
tioned in Section |Tj the need to reheat the Universe after 
inhation requires somehow a coupling between the SM 
and the inhationary sector, though there is no particu¬ 
lar constraint on this. Therefore, effective operators are 
expected to connect the Higgs with the inhaton when in¬ 
tegrating out some possible mediator held(s). However, 
as we will show in the following sections, the Higgs de¬ 
cays very fast after inhation into all SM species. Hence, 
even if there is an inhaton-Higgs effective coupling, we 
will assume in practice that its effect is negligible, with 
the possible Higgs-inhaton interactions not affecting the 
Higgs dynamics during or after inhation. 

Under these circumstances, the Higgs amplitude dur¬ 
ing inhation ’performs’ a random walk at superhorizon 
scales, reaching very quickly, within few e-folds, the equi¬ 
librium distribution [ 20 ] 


Peqip) = Afexp 


27r^ 




The correlation length, i.e. the physical scale above which 
the Higgs amplitude p huctuates according to Eq. (|^, 


^ Here we refer to a particle physics coupling, not the gravitational 
coupling. 


is given by ~ exp{3 .S/'s/a} [ 10], so it is ex¬ 
ponentially larger than the inhationary Hubble radius 
After the equilibrium distribution is reached at 
some point during inhation, the correlation length re¬ 
mains invariant until the end of the exponential expan¬ 
sion. Hence, immediately after inhation, the Higgs am¬ 
plitude p can be safely considered homogeneous within 
any volume of size I l^,. The Higgs amplitude varies 
randomly according to Eq. ([2), but only when we com¬ 
pare it at scales Z A, much larger than the correlation 
length. For convenience, we dehne 

a = X^p/H^f , (3) 

so that the distribution probability expressed over this 
dimensionless variable reads Peq oc exp[—27r^a"‘/3]. The 
roots of the moments of Peq are then given by 

c^ = {a^f^ = X^/\{p/HXf^, (4) 

where (...) denotes statistical average over the equilib¬ 
rium distribution in Eq. <§■ One hnds 

C 2 ~ 0.363, C 4 ~ 0.442, cg ~ 0.497, .... 

whereas 01=03 = 05 = ... = 0. We hnd a S [0.001,1] 
with 99.8% probability, whereas a < 0.001 holds only 
with a 0.17% probability, and a > 1 is yet further sup¬ 
pressed with a 0.03% probability. 

A typical amplitude of the Higgs at the end of inflation 
is given by its root mean square (rms) value 

¥>rms = C2P*/Al/4~1.15P,/Aj/t, (5) 

where we have defined the self-coupling normalized to a 
canonical value Ac = 0 . 01 , 

''^001 = = lOOA . ( 6 ) 

As we explain later, reasonable values of A are taken 
within the interval [10“^,10“®]. Hence, for A = 
10-2,10-3,10-4,10-5 (Ag^i = 1,0.562,0.326,0.178), we 
conclude that the typical Higgs amplitudes are of the or¬ 
der </5rms ~ independently of the value of A. We 
do not know the actual value of p within the ’progen¬ 
itor’ patch from which our visible Universe grew up. 
Actually, we do not know the value of the Higgs con¬ 
densate within any patch, we just know that typically 
p/H^ ^ 0(0.01) — 0(1) for reasonable values of A. That 
means that just after inflation, within any patch of size 
Z < Z*, the Higgs has a nonzero amplitude that could be 
really large, almost as big as 77* depending on its realiza¬ 
tion. The most updated upper bound for the inflationary 
Hubble rate is 

77* < 77i”“^ ~ 8.4 X lO^^GeV , 

so the Higgs amplitude at the end of inflation could be 
ranging around \p\ < (lO^^ — lO^^) QgV x (77*/77i™‘^’'^). 
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In order to analyze the dynamics of the Higgs after 
inflation, it is necessary first to fix the postinflationary 
expansion rate. Since we do not specify the nature of 
the inflationary sector here, we can parametrize the scale 
factor after inflation like 

/ 1 2 
a(t) = a* 1 +- <*) , p =——— (7) 

\ P ) (l + 3u;) 

with a* being the scale factor at the initial time t = 
(i.e. at the end of inflation), and w being the equa¬ 
tion of state of the Universe characterizing the expan¬ 
sion rate of the period following inflation. For instance, 
if the inflationary sector is described by an inflaton with a 
quadratic potential, the Universe expands as in a matter- 
domination (MD) regime during the inflaton oscillations 
following the end of inflation, so w = 0 and p = 2. If it is 
an inflaton with a quartic potential, the Universe expands 
as in a radiation-domination (RD) regime, with w = 1/3 
and p = 1. Since we do not really specify the inflaton 
sector, we are also free to consider other possibilities, 
including more ‘exotic’ scenarios where the background 
energy density decays faster than relativistic degrees of 
freedom, i.e. with w > 1/3 and p < 1. The paradigmatic 
example of this is a kination-domination (KD) regime, 
with w = 1 and p = 1 / 2 , obtained when an abrupt drop 
of the inflaton potential takes place at the end of infla¬ 
tion, transferring all the energy into kinetic degrees of 
freedom [301131] . 


A. Higgs oscillations 


The amplitude of the Higgs after the end of inflation is 
nonzero, and given that the Higgs potential is symmetric, 
the Higgs condensate is ’forced’ to oscillate around its 
minimum at (p = 0. The larger the Higgs amplitude, the 
sooner the oscillations will start after the end of inflation. 
The EOM (equation of motion) of the Higgs just after 
inflation is 


(p + 2'H(p + , ( 8 ) 


where ' = d/dt, and H is the comoving Hubble rate, given 
by 


n{t) = 


a 

a 


ci^ H^ 

[1 - t*)] 




(9) 


We will consider the evolution of the Higgs in an arbi¬ 
trary patch, inside which its amplitude [randomly drawn 
from Eq. & can be regarded as homogeneous. The cor¬ 
relation length is exponentially bigger compared to the 
Hubble radius, 4 ~ so if we just 

follow the Higgs within a causal domain of initial size 
I ^ l/i7» I*, then we can drop the Laplacian term 

on the rhs of Eq. (|^. The only scale in the problem is 
therefore so it is convenient to define a dimension¬ 
less conformal time z = — t*). We can then write 


the scale factor as a{z) = a*(l +p ^z)p. Introducing the 
variable 


h{z) 


a (f 

Oiijc 


( 10 ) 


with being the initial amplitude of the Higgs, we can 
rewrite the Higgs EOM in a more convenient form as 


32h3 = -h, 

a 


o2 _ 

p = — = VXa , 


( 11 ) 


where ' = d/dz, and /3 characterizes the frequency of 
oscillations. The term on the rhs scales as a”/a ~ 
(a*/a)^/^, and hence it becomes irrelevant very soon, 
since it decays as a"/a ~ z~^^p 1 . 

The initial condition for the Higgs amplitude in the 
new variables is, by construction, 


h* = 1 . 


( 12 ) 


The initial condition for the derivative h), = dh^,/dz = 1-1- 
(/(t*)/(a*i7»(^*), taking into account that the Higgs was 
in slow roll during inflation [i.e. </(t*) = — Aa^(/?2/27J*], 
reads out 

(13) 


The initial velocity of the Higgs and the frequency of its 
oscillations (in the dimensionless variables) both depend, 
through /3, on the initial amplitude of the Higgs (/?*, and 
the actual value of A. Therefore, at different patches 
of the Universe (separated at distances larger than the 
correlation length I ^ Z*), the Higgs will start oscillating 
with different amplitudes, and the oscillation frequency 
will also be different, see Fig. 

Depending on the amplitude of /?, the Higgs will start 
oscillating around the minimum of its potential sooner 
or later. This can be clearly seen in Eq. 0 , where 
the effective squared frequency of the oscillations of h{z) 
scales as oc /3^. For the canonical value of A = Ac = 
0.01 (Aooi = 1), the probability for the Higgs to start 
oscillating immediately at the end of inflation (i.e. that 
/? > 1 ) is extremely suppressed as 10 “^®^%, being even 
smaller for A < Ac (Aqoi < !)■ 

Therefore, at the end of inflation, the Higgs has, within 
any arbitrary patch of size smaller than Z*, an initial ve¬ 
locity in slow roll and a nonzero amplitude as large as 
p/Hif ^ 0(0.01)— 0(1). This amplitude remains ’frozen’ 
for a finite time until the start of the oscillations. Look¬ 
ing at Eq. Q, and denoting as Zosdfi) the time at which 
oscillations start at each patch, we see that the condi¬ 
tion for the onset of oscillations is a{zosc)^/Xp{zasc) = 
’H(zosc)- For simplicity, we will set the initial value 
of the scale factor to unity a* = a(t*) = 1, so that 
"H* = z = H^{t — t*), and a(z) = (1 -|- z/p)^. We 
will also denote any quantity evaluated at Zqsc with the 
suffix osc, so for example Oosc = a(-Zosc)- R follows that 


aoscVXposc = a-oscHosc = H^/alic, from which we find 
1 


V^^OSC — 


(o-osc) 


1 +- 


f/losc = l. (14) 
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taking into account that in our convention, a{zM) ^ 1- 
The period, in units of z, is then found to be 


Zt = 


7.416 7?* 

\/\(p{zM)a{zM) 


7.416 
PHzm) ' 


(15) 


Let us note that the factor 7.416 is only exact for RD. For 
MD or KD, one expects a similar though somewhat dif¬ 
ferent number, simply due to the term a"/a in Eq. 0, 
which affects the very early stages of the Higgs dynamics 
(even if it decays very fast after the onset of oscillations). 

We have obtained fits for Zqsc, hose, h{zM) and Zt as 
a function of /3 and for each postinflationary expansion 
rate, characterized by the equation of state u. These fits 
will turn out to be useful later on. We find at the onset 
of oscillations 

hose = 0.98/3-50^ (16) 


On the other hand, we find the field amplitude at z = 
Zm, and the oscillation period (measured from z = zm 
onwards), as 


/i(zm) = Ahose , Zt = 3(1+™) ^ (18) 


FIG. 2: Evolution of the Higgs field for /3 = 10“^, 2.5 x 
10“^, 5.Ox 10“^, 7.5x 10“^ and 10“'^ (corresponding to the red 
solid, orange dotted, blue dotted-dashed, green long-dashed 
and purple short-dashed lines, respectively). The background 
is RD, so w = 1/3. Dashed vertical lines mark the time 
Zosc(/ 3) when the oscillation condition is attained, = 77, 

whereas continuous vertical lines mark the time zm{P) when 
the first maximum in the oscillations is reached, characterized 
by the condition h'[ zm) = 0. Top: Evolution of h{z). Lower: 
Evolution of the physical Higgs which initially is frozen 

until the oscillations start, and then decreases as oc 1/a after¬ 
wards, as it oscillates. Similar plots are obtained for MD and 
KD backgrounds, whereas for other values of /3 the scale in 
the horizontal axis changes quite significantly. 


For a given expansion rate (characterized by the 
postinflationary equation of state w), the period of os¬ 
cillations depends sensitively on /3, since the period is 
fixed when the oscillation condition = 77 is at¬ 

tained at the time Zqsci which is itself a function of /? 
and w. The time scale zm at which h{z) reaches its first 
maximum, characterized by h'{zM) = 0, also depends 
consequently on /3 and w. The period of oscillation can 
be easily obtained from the case of a field with quartic po¬ 
tential, initial amplitude (p^, zero initial velocity </)* = 0, 
and RD background. In conformal time, when the scale 
factor at the onset oscillations is set to unity, it is given 
by r = 7.416 /(a/At^*) [S2]- In our case, we just need to 
count the oscillations from the first maximum at z = zm, 


where {A,B) ~ (1.28,6.30), (1.22,6.25), (1.17,6.25) for 
w = 0, 1/3, and 1, respectively. 

At the end of inflation, the Higgs energy density at a 
given patch is mostly dominated by its potential energy, 

K = ^ (19) 

which represents a very small contribution of the total 
energy budget at that moment, p* = 3mpi7^. Averaging 
over realizations, we find 


(K) 4 

3^2772 C4 



~ 4 X 10”^^ 


77, 

^{max) 


2 


( 20 ) 


At the onset of oscillations, part of the potential energy 
will become kinetic, with the two contributions - kinetic 
and potential - becoming of the same magnitude. In 
order to see this, let us first write the total energy density 
of the Higgs as 


with the kinetic and potential contributions given by 

= EKiz,P) + Ev{z,p). (22) 
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FIG. 3: Top: Total energy p^jV, (continuous lines) and 
its oscillation-averaged value ptp/V* (dashed lines), for /3 = 
10“^2.5 X 10■^5.0 X 10■^7.5 x 10■^ 10“^ (from right 
to left, red, orange, blue, green, purple). Vertical grey 
lines mark zm{P), signaling from what point the averaged 
curves should be considered valid. Bottom: The func¬ 
tions Ek{z,P) (dashed, purple), Ev{z,fi) (dotted-dashed, 
blue) and E{z, j3) = Ek{z,I 3) + Eviz, j3) (solid, red), and 
their averages Ek{P) (purple), Ev{/3) (blue) and E{/3) = 
Ek{P) + Ev{P) (red), for /3 = 10“^. The figure shows that 
EviP) = \Ek{I3) = \E{P). All plots obtained for a RD 
background. 


We can then take the average over the Higgs oscillations 


(23) 

_ 1 rz+ZT{0) 

= (24) 

and again split the result into potential and kinetic con- 


^ Note that we are not inclu ding; in the average the prefactor 
l/a'^(t) factorized out in Eq. ( |2l[ |, since the scale factor changes 
only marginally during each oscillation. Therefore, we are only 
averaging the contribution due to the Higgs oscillatory behavior. 


tributions, i?(/3) = Ek{I 3) + Ev{l3), where 


EviP) 

EK{f3) 


1 rz+ZrriP) 1 _ 

ZW) L 


a',Y 2^, 


1 PZ + ZTiP) _ 

(26) 


In Fig. we can see the total energy density of the Higgs 
for different values of /3, with and without averaging. Of 
course, the oscillation-averaged expressions are only valid 
when the Higgs has started oscillating at z > zm, as 
clearly appreciated in the plot. The figure also shows 
very nicely the fact that the averaged components verify 
Ev(P) = \E{I3) and Ek{P) = ^E{f3). Possibly, the 
most relevant aspect to be remarked is the well-known 
fact that the Higgs energy density scales as with the 
expansion of the Universe |33] . behaving as if it were a 
fluid of relativistic species. 


III. HIGGS DECAY: ANALYTICAL ESTIMATES 

As just explained, the Higgs oscillates everywhere in 
the Universe, although the time to start the oscillations 
depends sensitively on the initial condensate amplitude, 
which varies from patch to patch according to Peq{f)- 
Once the oscillations have begun within a given patch, 
all fields coupled directly to the Higgs are excited every 
time the Higgs goes through the minimum of its poten¬ 
tial. In the case of bosonic species, this is known as para¬ 
metric resonance, since a cumulative effect takes place, 
producing a resonant growth of the number density of 
species O [TOl |32l [34U37] . Although there is no paramet¬ 
ric resonance in the case of fermionic species, yet an inter¬ 
esting effect occurs, since modes with successively higher 
momenta are excited as the oscillations carry on [55VH5] . 
For a review of parametric excitation of fields in the sim¬ 
ilar context of preheating, see [HI SS] . 

All charged leptons of the SM are directly coupled 
to the Higgs via a Yukawa interaction, so all fermions 
of the SM will be excited during the oscillations of 
the Higgs [43], with the possible exception of neutrinos. 
Among the SM fermions, the top quark has the largest 
coupling to the Higgs, so most of the energy transferred 
into fermions goes into top quarks. More importantly, 
the SU{2)l gauge bosons are also coupled directly to 
the Higgs, and indeed the strength of their coupling is 
very similar to that of the Yukawa top quark. When two 
species, one fermionic and another bosonic, are coupled 
with the same strength to an oscillatory homogeneous 
field, the first burst of particle production is actually spin 
independent, and hence an equal number of bosons and 
fermions are created [46j. However, the successive parti¬ 
cle creation bursts at each Higgs zero crossing take place 
on top of an already existing number density of previ¬ 
ously created species. The spin statistics becomes then 
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crucial, differentiating bosons from fermions in a notice¬ 
able way: bosonic occupation numbers start growing ex¬ 
ponentially as the oscillations accumulate, whereas the 
fermion occupation numbers are always Pauli-blocked, 
forcing the transfer of energy into modes with higher and 
higher momenta. Both bosonic and fermionic excitations 
represent a sizable transfer of energy from the Higgs con¬ 
densate. However, for equal coupling strength [as it is the 
case between top quarks and SU{2)l gauge bosons], the 
transfer of energy is much more efficient into the bosonic 
species [ID]. Besides, in the context under study here - 
the decay of the Higgs after inflation -, the subdominant 
production of the SM charged leptons has been already 
addressed in [43] . Therefore, in this paper we will only 
focus on the production of the most energetically dom¬ 
inant species among the Higgs decay products, the 
and Z gauge bosons. 

In order to study the dynamics after inflation of the 
Higgs and its most energetic decay products, one should 
in principle consider the full SU{2) x 17(1) gauge struc¬ 
ture of the SM electroweak sector. However, one can 
make reasonable approximations for both analytical and 
computational purposes. In this work we have consid¬ 
ered the following approximate schemes, mimicking the 
structure of the SM interactions: 


The approach i is our most precise modeling of the SM 
interactions, but also the most involved one. We thus 
postpone its implementation for later on in Section jV] 
The approach ii, though less accurate, has a clear advan¬ 
tage versus the gauge case: it allows not only for a lattice 
implementation (which we introduce in Section IV), but 
also for an analytical treatment (which we present in the 
remaining of this section). The analytical estimates rep¬ 
resent only an approximation to the system described by 
the scenario ii, but yet provide a valuable insight into 
the understanding of the dynamics. The order of pre¬ 
sentation in the paper of our different approaches is thus 
based on increasing progressively the degree of proxim¬ 
ity to the real system. First, in the remainder of this 
section, we start with the analytical treatment of the 
global modeling, ignoring all nonlinearities of the sys¬ 
tem. In Section jlVj we implement the global scenario 
ii) on the lattice. This way, we fully capture all nonlin¬ 
earities within this modeling, even if we yet neglect the 
gauge nature of the interactions. Finally, in Section jV] 
we present a lattice implementation of an Abelian mod¬ 
eling of the system. This fully captures the nonlinearities 
within such modeling, while preserving at the same time 
the gauge-invariant nature of the interactions. 


i) Abelian model. This consists in modeling the in¬ 
teractions between the electroweak gauge bosons 
and the Higgs with an Abelian-Higgs analogue. 
Since gauge fields are initially excited by the Higgs 
from the vacuum, it is clear that nonlinearities due 
to the truly non-Abelian nature of SU{2) are ex¬ 
pected to be negligible during the initial growth 
of the gauge field occupation numbers [17]. The 
authors of Ref. m have shown that using the 
Hartree approximation, the effective contribution 
induced by the created gauge bosons onto them¬ 
selves (due to the non-Abelian nonlinearities) can 
be neglected as long as the backreaction from the 
gauge fields onto the Higgs does not become signif¬ 
icant. In principle, this fact fully justifies ignoring 
the non-Abelian structure of the SM interactions, 
while maintaining only the Abelian dominant part. 

ii) Global model. A more crude approximation can 
yet be done, by ignoring the gauge structure of 
the interactions. This does not mean that we ig¬ 
nore the interactions themselves, but rather that we 
consider them as if they were dictated by a global 
symmetry, instead of a gauge one. In this scenario, 
one simply solves the mode equations of various 
scalar fields coupled to the Higgs with a quadratic 
interaction. Each of these scalar fields mimics a 
component of the gauge fields, with the quadratic 
interactions reproducing the coupling of the gauge 
bosons and the Higgs obtained from the SM gauge 
covariant derivative terms. This way, one can pre¬ 
sumably capture the initial stages of the parametric 
resonance of and Z bosons. 


A. Analytical approach to the Higgs decay 

In principle, we can follow the initial stages of the para¬ 
metric resonance of the , Z bosons by simply solving 
the mode equation for a scalar field X) coupled to the 
Higgs with an interaction term of the form -^ci- 

alytical results following this approach have indeed been 
presented in |D] , so our work in this section should be un¬ 
derstood only as complementary to such reference. We 
develop nevertheless some new formulas which will be 
useful later on, in order to assess the reliability of this an¬ 
alytical approximation when compared to the fully non¬ 
linear numerical lattice simulations. 

The equation for the Fourier modes of the field Xj after 
an appropriate conformal redefinition Xfe = Af^/a, and 
assuming RD, can be mapped into [32] 

A" + (k2 + q{h/hosc)^) Afc = 0 , 9 = y > (27) 

with q being the resonance parameter, k = k/{'/Xiposc), 
' = d/dz, and z = 'Hosct- Given the behavior of h(z), dic¬ 
tated by the Higgs quartic potential, this equation cor¬ 
responds indeed to the Lame equation [32] . which has 
a well-understood structure of resonances. Whenever 
q S 7[n(n-|-I), (n-|-I)(n-|-2)], with n = 1, 3, 5,... (i.e. q G 
[I, 3], [6, 10], ...), there is an infrared band of resonance 
0 < fc < fc* = for which Xk cc with Ufc 

bigger the smaller the k (therefore maximum at k = 0). 
If the resonance parameter q > 1 is not within one of the 
resonant bands, but lies in between two adjacent bands, 
then there is still a resonance of the type Xk oc but 



within a shorter range of momenta 0 < fcmin ^ k 
and hence with a smaller Floquet index There is a 
theoretical maximum value for the Floquet index given 
by = 0.2377... [35], so that any fik is always con¬ 
strained as /ifc < for g > 1. For resonant parame¬ 

ters g 3> 1, is typically of order ^ 0(0.1); see Fig. 

For simplicity, we will consider until the end of this 
section that the resonance parameter g = e^/A always 
falls within one of the resonant bands, g G [1,3], [6,10], 
[15,21],.... As a matter of fact, in order to identify 
with the gauge coupling between the Higgs and a 
gauge field, we need to make the identification —>■ g^/4, 

with the gauge coupling g| or of either the Z or 
the gauge bosons. This matches correctly the in¬ 
teraction derived from the covariant gauge derivative of 
the electroweak sector of the SM. The gauge couplings 
of the Z and gauge bosons verify g| « 2g^ « 0.6 
at very high energies. Due to this relation, it is likely 
that either qw = g^jiX or qz = gz/^X « 2gw, will fall 
within one of the instability bands. Let us note, how¬ 
ever, that we cannot predict this, since the value of the 
Higgs self-coupling A at high energies is quite sensitive to 
the uncertainties in the Higgs mass mu-, the top quark 
mass mt, and the strong coupling constant as- Conse¬ 
quently, we cannot really know the exact value of these 
resonance parameters. However, in order to guarantee 
that during inflation the Higgs fluctuations remain be¬ 
low the critical scale /r_|_ (above which the self-coupling 
starts decreasing, dX/dg, < 0), and taking into account 
that the inflationary Hubble rate is constrained from 
above as i7* < 10^"^ GeV [29], the Higgs self-coupling 
value at high energies can then only be within the range 
10“^ ^ A < 10“®. Pushing A to smaller values is in prin¬ 
ciple possible, but it represents a fine-tuning and requires 
some of the parameters mt,mH, 0 's to be more than 3 
sigma away from their central values. We will consider 
therefore the range 10“^ ^ A < 10“® as the only ac¬ 
ceptable one (with A ^ 10“® only marginally valid). If 
beyond the SM physics affects the running, say stabiliz¬ 
ing the Higgs potential at high energies, then A remains 
positive and typically of the order A ~ 10“^ —10“^. Con¬ 
sidering the range 10“^ ^ A < 10“®, and taking into ac¬ 
count the strength of the , .Z gauge couplings at high 
energies, we obtain that the resonant parameters can 
only possibly be within the range 0(10) < g < O(IO^). 
In particular, since at high energies g^ = ~ 0.3 

for W gauge bosons, we obtain g = 7.5 for A = 10“^, 
and g = 3000 for A = 2.5 x 10“^. For Z bosons we 
obtain similar resonance parameters, but twice as big. 
For completeness, we have sampled resonance parame¬ 
ters within the interval g G [5,3000], which corresponds 
to a range A = 1.5 x 10“^ — 2.5 x 10“^ for W bosons and 
A = 3.0 X 10“^ — 5.0 X 10“® for Z bosons. 

Let us then consider just one particle species rep¬ 
resenting either the Z or one of the W gauge bosons, 
that will be parametrically excited during the Higgs os¬ 
cillations. Let g^ be the coupling strength to the Higgs 
and let us represent the gauge field as if it were simply a 


collection of three scalar fields (one for each spatial com¬ 
ponent), all coupled with the same strength to the Higgs. 
The growth of the fluctuations in the initial stages of res¬ 
onance is described by the linearized Eq. (27). As long 


as the linear regime holds, even if the amplitude of the 
fluctuations grows exponentially, the use of three scalars 
should represent a good mapping of the real problem of 
gauge field excitation. Of course, one is ignoring this way 
the backreaction of the created bosons into the Higgs, as 
well as certain contributions in the gauge fields’ EOM, 
which should be present if the gauge symmetry was re¬ 
stored. 

The energy density of the created particles due to the 
resonance is then given by 


PA = 


27r2o^ 


r ^2 2 _ 

/ dkk'^UkUJk , ujI = -^ + ^ ( 28 ) 


with the factor 3 accounting for the three spatial compo¬ 
nents of a gauge field, and where we have introduced an 
oscillation-averaged effective mass for the gauge boson. 


2 


gW. 


I 


gZ + ZTi0) 


4 a2 Zt{/3) 


dz'h?{z') . (29) 


For g 3> 1, the maximum (comoving) momentum possi¬ 
bly excited in broad resonance is given by 

gl/2 „i/2 

kl = , (30) 

from which, given that ^ we see that 

-^~O(10)gi/2»l. (31) 

(k^/ay 

In other words, in broad resonance g 1, the de¬ 
cay products are always nonrelativistic, and correspond¬ 
ingly we can approximate the effective mode frequency 

as Wfe ~ niA ~ where /irms = V^. It 

turns out that /irms — hose independently^ of j3. If g 
is within a resonant band, then all modes with momenta 
G y k < kit are excited with some Floquet index vary¬ 
ing within [0, g^“’^’'^(g)]. This corresponds to the cases 
with blue solid lines in Fig. [^ We can therefore model 
the occupation number of the excited modes simply as a 
step function Uk = e2Ffe!^0(l — A:/fc*), where ~ 0(0.1) 
and y = hLoscid ^osc) — aoscV^A(^osc/77.,c)(z 2^osc) — 


® For /3 > 0.3 there is some dependence, but still '/h^/hosc ~ 
0 ( 1 ). 
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where we have used the fact that /3/iosc = l/^Oosc- 

This is how the energy density of the gauge bosons 
(those fully within a resonant band) will grow, at least as 
long as their backreaction into themselves and/or into the 
Higgs remains negligible. Using this linear approxima¬ 
tion we can estimate the moment ZeS at which an efficient 
transfer of energy has taken place from the Higgs into the 
gauge bosons, characterized by pAizes) = Pipizes)- This 
will be just a crude estimate of the time scale of the 
Higgs decay, since by then backreaction and rescatter¬ 
ing effects have become important, invalidating the linear 
approach. However, the nonlinear effects due to backre- 
action/rescattering of the decay products, simply tend 
to shut off the resonance. Hence, using the linear regime 
for inferring the Higgs time scale should provide, at least, 
a reasonable estimate of the order of magnitude. More 
importantly, it provides the parametric dependences of 
both the time when the resonance is switched off, and 
the moment when the energy has been efficiently trans¬ 
ferred into the gauge bosons. 

The energy of the Higgs, since the onset of the oscilla¬ 
tions, decays as 


= ^{h/h,sc) \ > (33) 

4A (a^/On 


where (h/hoscY ^ 0{1). We can now find by simply 
equating Eqs. (32) and (33), 


U/4 


25/27r5 

so that 


{h/hoscYe 


2p = ^(^h/hoscY, (34) 


Y ^osc 

— -^^osc “r ^ _ 


log 


jh/hoscY 

iYrras/h 

osc) 


(35) 


k/ (VTt/Josc) 



3- 25/27r5\ 
5^ / 



FIG. 4: This shows the band structures of the Lame equation, 
Eq. (271, for several resonance parameters ranging between 
q = 5 and q = 3000. In each panel, we plot the corresponding 
Floquet index (where oc as a function of the 

momentum k. We can divide the different q into two groups: 
those which contain a resonance at ft = 0 (blue lines) and 
those which not (purple dashed lines). 


Let us recall that ~ 0.3,0.6 at large energies, and 
q = g^/(4A) ^ 0(10) — O(IO^), depending on the value 
of A. Taking this into account, we find that the first term 
in the brackets of the rhs is always irrelevant, the second 
term is constant and of the order ~ 9, and the last term 
is of order ~ —1. Therefore, we can approximate the 
above expression, using = (1 -I- ^Zqsc), as 


(flosc) 'p {z — Zosc)) with z = H^t. It follows that 



(32) 


— -^osc “t“ ^_ 


Pk 


It _ I 2:osc 

PPk. 


PPk 


(36) 


Looking at Fig. we see that the Floquet index of the 
modes 0 < fc T fc* for which q is within a resonant band 
(blue solid lines of the figure), can be well approximated 
by a simple step function pk ~ Pk^Y ~ k/kY), with a 
























































































10 


mean Floquet index Jlf, ~ 0.2. Taking this into account 
and using the fit of Eq. (16) for the time scale at the 
onset of oscillations ZoscipYT^e hud 


-^eff 


20 X 



/?■ 


(1 + 3^.) 
3(1 + ™) . 


(37) 


The scale factor at z = Zeff is then given by 

Oeff = a(zeff) (20(1 + 3w)) (1+3™) . /3“ 3(1+™) _ (|38) 


It is clear that depending on how small the initial 
value of P is within a given path of the Universe, the 
longer it takes for the Higgs to transfer energy efficiently 
into the gauge bosons, simply because the longer it takes 
(since the end of inflation) to start oscillating. Since 
/3rms 0(0.1), we see that typically the Higgs decays 
at a time Zeff(/3rms) ~ O(IO^). Although the time varies 
from patch to patch depending on the values of /3, it 
is clear that the Higgs tends to decay really fast after 
inflation, within a few dozens of oscillations. In the fol¬ 
lowing sections we will check the validity of this estimate 
by comparing it with the outcome obtained directly from 
lattice simulations. 


IV. LATTICE SIMULATIONS, PART 1: 
GLOBAL MODELING 


In this section, we continue modeling the SM inter¬ 
actions with a set of scalar fields. More specifically, we 
consider the Lagrangian 




xl, 


(39) 


with i = 1,2, 3. Varying the action S = J dPx C leads to 
the classical EOM 

0 + 2 %^- VV + a"(AX + e' E «> (40) 

i 

Xt + 2'Hxi - = 0 • (41) 


The term e^^p^Xi 


under the identification = g^/4, 
mimics precisely the interaction term from the covariant 
derivative of the EW gauge bosons, where 

stands for either or W^, and $ is the Higgs dou¬ 
blet. More concretely, choosing the unitary gauge for the 
Higgs $ = (0, and fixing Aq — 0, we can identify 

each Xi with each spatial component of the gauge boson 
Ai, and p with the unitary representation of the Higgs. 
This way, by solving the system of scalar held equations 
(40) and (41), we can study the properties of the Higgs 


interactions with gauge bosons in an approximative way. 

In Section pH A| we studied this scenario, following the 
huctuations of the helds Xi with the help of the analytical 
solutions of the Lame equation, Eq. ( |^ . We exploited 
the band structure of this equation and used some ap¬ 
proximations in order to arrive at our analytical results. 


summarized in Eqs. (35 )-(|38[). In reality, the scalar helds 


Xi follow the Lame equation only initially, in the regime 
when the nonlinearities (due to their small backreaction 
onto the Higgs) can be neglected. The huctuations of 
the Xi helds grow exponentially during the linear regime, 
and as we will show, it does not take long until they 
start to impact onto the Higgs dynamics. At that mo¬ 
ment, the system becomes nonlinear, and only by fol¬ 
lowing in parallel the coupled EOM of the Higgs and Xi 
helds, can we really understand the held dynamics within 
this modeling. The aim of this section is, therefore, to 
solve numerically in a three-dimensional lattice the sys¬ 
tem of equations (|40|) and (41). Only in that way can 


we fully capture the nonlinear behavior of this system 
beyond the simpler linear regime of the Lame equation. 

We now present the main results of the lattice simu¬ 
lations carried out for this scenario. We start with the 
following change of held variables: 




y _ Xi a 

WTT ■ 

-n. * ct* 


(42) 


It is also convenient to redehne new spacetime coordi¬ 
nates z^ = (z°,z®) with respect to the conformal ones 
= (a:°,a;*) = {t,x''), as 

z = z^ = HA , 2 * = . (43) 

With these held and coordinate redehnitions, we elimi¬ 


nate the friction terms in Eqs. (40) and (41), and pro¬ 


duce an equivalent set of dynamical equations, written in 
terms of the new dimensionless variables: 


h” - + e^h 


j 


X'l - + qp^ypx, = 


“A 
= —h, 

a 

(44) 

“"a 

(45) 


with ' = d/dz, and the spatial derivatives taken with 
respect to the z* variables. A lattice version of these 
equations is presented in Appendix As already men¬ 
tioned, we will identify —)■ g^/4, with being either 
or g^. The resonance parameter that appears natu- 
J L ^ 

rally in Eq. (45), q = ^, should therefore be interpreted 


as(Z= fy. 


We have solved Eqs. (44) and (45) in three-dimensional 


lattices with periodic boundary conditions. We consider 
initial conditions given by a homogeneous Higgs mode 
(as described in Section 


h(o) = i, hioy = i-^. 


(46) 


and a null zero mode for the scalar helds coupled to the 
Higgs, 


A,(0)=0, A'(0)=0. 


(47) 


We add, on top of the homogeneous contributions, a set 
of Fourier modes with spectrum (|//cp) = 2 a?uik (™ p4^y®“ 
ical variables), mimicking the quantum vacuum huctua¬ 
tions of the ground state of a scalar held in a FRW back¬ 
ground. Let us recall that the Higgs is frozen in slow 
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roll until the oscillation condition Eq. (14) is attained at 
■z = 2 osc; see the bottom panel of Fig. Hence, dur¬ 
ing the time 0 < z < Zqso "we only evolve in the lattice 
Eq. (44), corresponding to the slow rolling of the Higgs 
field (tire homogeneous mode of the Xi fields is kept to 
zero). At z = Zosc, we add the small inhomogeneous 
Fourier modes to all fields, and from then on, we evolve 
together Eqs. ( [44| ) and ( [4^ . The reader can find more 
details about our methodology for introducing the initial 
conditions in Appendix [B| 

We have run simulations for different values of /3 = 
Since we know that 10 ^ < a < 1 and 
10“® < A < 10“^, we find that 10“^ ^ /3 < 1- We 
have thus run simulations for /? = 0.5, 0.1,10“^, 10“^ 
and 10“"^. Note that the root mean square of (3 is 
^ 0.115Ao^i « 0.115,0.065,0.037,0.020 for Aqoi = 
1,10“^, 10“^, 10“^, respectively. The probability distri¬ 
bution Eq. <111 for a (and hence /?) is very non-Gaussian 
and, independently of A, /3 > 0.5 is exponentially sup¬ 
pressed. The range 10“^ < j3 < 0.1 is obtained with 
more than 99% (99.8%, 99.7%, 99.4%, 99.1% for Aooi = 
1 , 0 . 1 , 0 . 01 , 0 . 001 ), while /3 < 10 “"^ is attained with < 1 % 
for all values of A. Hence the values of /? that we have 
chosen, /3 G [10“"*, 0.5], sample fairly the range of random 
initial Higgs amplitudes dictated by Peq( 7 ’)- 

The actual value of A is quite uncertain, since it de¬ 
pends on the energy scale of inflation. Besides, for a 
given Hubble rate 77*, it can still vary significantly given 
the uncertainties in mff,as and mt (mostly in the lat¬ 
ter). Due to this, for each value of we have chosen a 

2 

set of 26 resonance parameters q = ^, logarithmically 
spaced between q = 5 and g = 3000. This corresponds 
to sampling the Higgs self-coupling from A ^ 10“® to 
A ^ 10“^. Scanning this way /3 and q led us to charac¬ 
terize the behavior of the system, scrutinizing all possible 
different outcomes depending on A and </?*. In Table [ij we 
list the values of all the resonance parameters q that we 
have considered. We have guaranteed that by sampling 
different values, we include both the cases in which q is 
within a resonance band of the Lame equation, or in the 
middle of two bands (see Section HI). 

Note that we have run simulations for three different 
expansion rates, corresponding to a MD universe, a RD 
universe, and a KD universe, given by w = 0, | and 1 in 
Eq. 0, respectively. The following results in this section 
will be presented only for a RD background. The gen¬ 
eralization to other expansion rates will be considered in 
Section IVII 

Our simulations depend only on two parameters, q and 
/I. For each pair of values (q, /3), we have run simulations 
on a lattice with A7 = 128 points per dimension, with 
periodic boundary conditions. The minimum momentum 
captured in each run is km = with dx being the 

lattice spacing. The maximum momentum sampled in 
the lattice is The length of the lattice box 

side is L = Ndx. For each value of /3 and q, we have 
made sure that our results are not sensitive to the lattice 
spacing dx and/or the lattice size L. More details about 


qw 

Aooi 

^min i^QW ) 

^max(QW ) 

5 

1.5 

0.72 

1.09 

6 

1.25 

0 

0.97 

8 

0.938 

0 

0.69 

9 

0.833 

0 

0.49 

11 

0.681 

1.33 

1.54 

14 

0.536 

0.67 

1.26 

18 

0.417 

0 

0.83 

23 

0.326 

1.43 

1.72 

29 

0.259 

0 

1.24 

37 

0.203 

1.75 

2.02 

48 

0.156 

0 

1.22 

61 

0.123 

1.36 

1.92 

79 

0.095 

2.06 

2.38 

101 

0.074 

0 

0.91 

130 

0.058 

0 

1.10 

167 

0.045 

0 

0.88 

214 

0.035 

2.34 

2.83 

275 

0.027 

0.56 

2.21 

354 

0.021 

2.71 

3.22 

454 

0.017 

0 

1.43 

584 

0.013 

0 

1.42 

750 

0.010 

2.93 

3.65 

1030 

0.0073 

1.18 

3.04 

1550 

0.0048 

0 

2.85 

2200 

0.0034 

0 

1.37 

3000 

0.0025 

0.87 

3.72 


TABLE I: Different resonance parameters q used in the sim¬ 
ulations, together with the corresponding values of the Higgs 
self-coupling derived for ~ 0.3. For each case, we 

also provide the minimum and maximum momenta (in units 
of H,), kmin < k < fcmax, of the first resonance band. Half 
of the cases have a band down to fcmin = 0, while the others 

have fcniin > 0. 


these issues are given in Appendix |A| 

In Fig. we plot, as a function of time, the volume- 
average of the modulus of the (conformally transformed) 
Higgs field \h\. In this figure, we show the outcome cor¬ 
responding to /3 = 0.01, and four different resonance pa¬ 
rameters, 9 = 8, 14, 101 and 354. The values q = 8 ,101 
are centered close to the middle of a resonance band of 
the Lame equation, while q = 14, 354 are between adja¬ 
cent bands. In this figure we also show the corresponding 
envelope curve of the Higgs oscillations. One conclusion 
is immediately clear: the time scale of the Higgs ampli¬ 
tude decay depends noticeably on q. By running simula¬ 
tions for each of the q values displayed in Table [I) we have 
fully characterized the q dependence of the Higgs decay. 
Note that in Table |T] we have also indicated the range 
of momenta fcmin 33 k < fcmax excited for each value of 
9 , according to the Lame equation. Such a range corre¬ 
sponds to the band with the largest Floquet index p.niaxj 
which coincides in all cases with the most infrared band; 
see Fig. The ^{k) index was obtained by solving the 
Lame equation for a given q parameter, and finding the 
range of momenta such that /r(fc) > 0. The band struc¬ 
ture can be well appreciated in Fig.|^ where we plot ^(fc) 
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FIG. 5: Volume-averaged value of the Higgs field \h\ as a 
function of time, for four different resonance parameters, q = 
8,14,101 and 354. Also plotted, the corresponding envelope 
functions of the oscillations. All cases correspond to j3 = 0.01. 


for each of the values of q listed in Table I. As mentioned, 
we have sampled all possible cases, including when q is 
within a resonant band (either close to the middle or to 
the extremes of the band), and hence fcmin = 0, or simply 
outside of any band (between adjacent bands), and then 

^min ^ 0- 


Before examining in more detail the general behavior 
of all the fields in the system, we can make some com¬ 
ments about the Higgs behavior. First of all, let us note 
that h oscillates with a period T which is, as expected, 
independent of the value of q. Even if it cannot be really 
appreciated in Fig.[^ we have checked that the period co¬ 
incides initially with the analytical expression given by 
Eq. (15), until it becomes slightly modulated due to the 
interactions with the y fields (though it does not change 
significantly). Looking at the different panels of Fig. it 
seems that the Higgs decay is slower the greater the res¬ 
onance parameter q is. This is very opposite to the intu¬ 
ition gained by the study of the Lame equation in Section 
|HH which dictates that the larger the q, the shorter the 
decay time of the Higgs.^ We thus see on this the first 
difference between the simplified study of the system of 
scalar fields in the linear regime (Section [III A ), and the 
real outcome when nonlinearities are incorporated in lat¬ 
tice simulations. We will further comment on this issue 
later on. 


One can distinguish two different stages in each decay 
process. Let us look, for instance, at the upper panel of 
Fig.ll where the Higgs modulus is plotted for q — 23, 


^ Contrary to ’popular wisdom’ about parametric resonance, the 
time scale identified with the ’oscillatory field’ decay time 
in the linear analytical approximation, is in practice mostly in¬ 
dependent of q. It is true that the larger the q the shorter the 
decay, but the dependence is only logarithmic [recall Eq. l |35[ l], 
and the number of oscillations does not change appreciably. 




FIG. 6: Volume-averaged value of the Higgs modulus for q = 
23, P = 0.01 and RD. An initial plateau until z = Zi can 
be clearly distinguished in the top panel, where we plot the 
conformally transformed Higgs. At later times 2 ; > Zi, the 
amplitude of the Higgs drops abruptly, due to its decay into 
the X fields. In the lower panel we plot the physical Higgs 
= \h\/a, where we can appreciate that the plateau for 
h translates into a dilution oc 1/a for ip, due to the expansion 
of the universe. The decay of the Higgs into the other fields 
at later times is manifested by a significant decrement of |yj| 
well below the 1/a decaying envelope. 


and where we also include the envelope curve of the oscil¬ 
lations. One can clearly appreciate that initially, and for 
some time, the envelope is approximately constant, re¬ 
ducing its amplitude only slightly. This is observed as a 
plateau feature in the upper panel of Fig. The vertical 
dashed line in the figure indicates the end of this ini¬ 
tial behavior, after which a second stage of rapid decay 
follows. Let us note that when we talk about the de¬ 
cay of the Higgs amplitude, we refer to the conformally 
transformed one h. The amplitude of the physical Higgs 
p/p* = h/a{t) is always decaying with the scale factor, 
no matter what its coupling to other species is. Before 
the second stage starts, the physical Higgs amplitude p 
decays mostly due to the expansion of the Universe, and 
not because of an efficient transfer of energy into the 
scalars. However, both effects are combined afterwards, 
producing an even more sharp decay of the physical am¬ 
plitude. This is clearly seen in the lower panel of Fig. 
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In order to understand better this two-stage behavior, 
we plot the different contributions to the total energy of 
the system as a function of time. The energy density can 
be conveniently written as 


p{z) = 14 


Et{z) 
a{zY ’ 



Et{z) — E‘^ + Ey + Eq + E^ -|- Eq + i?int , 


(48) 

(49) 


where, for our choice of variables, the Higgs and x field 
contributions to the kinetic (K) energy are given by (' = 
d/dt, ' = d/dz) 



^ 14 2a2 13^ \ a 


jex - a^X^X^ 

K - K 2a2 


2A 




the gradient (G) contributions by 

,2 


O - u 2 a 2 


K 2 a 2 




(50) 

(51) 


(52) 

(53) 


i=l 



and finally, the Higgs potential (V) energy and the inter¬ 
action (int) term, by 


Ey = 
E'lnt = 


14 4 




2 e 2 


Y 2 X*Xi ~ ^2^ ■ 


(54) 

(55) 


In Fig. we have plotted the different contributions 
to Et{z) for the parameters /I = 0.01 and q = 8. Ini¬ 
tially, the system is dominated by the kinetic and poten¬ 
tial energy densities of the Higgs. This corresponds to 
the regime of anharmonic oscillations of the Higgs con¬ 
densate described in Section |TTj for when the coupling to 
other fields was ignored {g^ —> 0). However, in reality, 
as soon as the Higgs starts to oscillate, there is an en¬ 
ergy transfer into any species coupled to the Higgs. Each 
time the Higgs crosses zero, a fraction of its energy goes 
into the x fields. Initially, the amount of energy trans¬ 
ferred at each zero crossing is small relative to the to¬ 
tal energy stored in the Higgs. Therefore, it takes some 
time until the transfer becomes noticeable. The Higgs 
energy components represent the dominant contribution 
to the total energy during the initial oscillations, so the 
Higgs evolves initially without really noticing the pres¬ 
ence of the other fields. Eventually, at the time z = zi, 
the energy transferred into the x fields becomes signif¬ 
icant enough, compared to the Higgs energy itself, say 
a fraction p^jPip = 5 < 1. The Higgs condensate be¬ 
comes affected by the transfer of energy into the other 
fields when 8 > S{zi) = 0.1. From then onwards, the 


FIG. 7: Top: We show the envelope curves of the oscilla¬ 
tions of the different contributions to the total energy Et{z), 
obtained for q — 8, jd = 0.01 and RD. The gray, vertical 
dashed line corresponds to the decay time for q = 8. Bot¬ 
tom: Same quantities as in the upper figure (same color cod¬ 
ing), but zooming in the area of interest. We also add two 
new lines, a pink one corresponding to the sum of the Higgs 
gradient energy and the interaction energy, and a light blue 
line, representing the sum of the y fields’ gradient energy plus 
the interaction energy. We see that the decay time indicates 
equally good the time when the Higgs kinetic energy stops 
decaying, and the time when equipartition is set. 


Higgs continues pumping energy into the other fields at 
z > Zi, but the amount of energy transferred at each zero 
crossing is no longer a small fraction of the energy avail¬ 
able in the Higgs condensate itself. Therefore, soon after 
backreaction becomes noticeable at z = Zi, the previ¬ 
ously exponential growth of the y fields energy densities 
stops, eventually saturating to a fixed amplitude. This 
is clearly seen in Fig. [Tj where the gradient and kinetic 
energy densities of the y fields saturate to an almost con¬ 
stant amplitude. This happens because the Higgs has not 
enough energy anymore to accomplish transferring a siz¬ 
able fraction of energy into the y fields. At the same 
time, immediately after z = Zi, the Higgs energy density 
drops abruptly. This is so because the amount of energy 
transferred from the Higgs into the other fields, even if 
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FIG. 8: Top: The different times Zi{q) obtained from our 
simulations, for RD {uj = 1/3) and /3 = 0.01. Purple triangles 
and blue circles correspond to q parameters inside or outside 
a resonance band of the Lame equation respectively. The blue 
and purple continuous lines correspond to the best fit to the 
circles and triangles respectively; Eq. (561. The dashed line 
corresponds to the analytical estimate z^g ~ 200, obtained 
from Eq. (371 with /ife = 0.2. Bottom: The different points 
show the Higgs time decay Ze{q) as a function of q obtained 
from our simulations for the same (a;, /3) values as the upper 
panel. The brown line corresponds to the best fit, Eq. (571. 


not significant anymore compared to the energy stored 
in the x fields (hence the saturation of their growth), 
represents a significant fraction of the energy available 
in the Higgs at that moment. Therefore, the energy of 
the Higgs (mostly dominated by the kinetic contribution) 
drops abruptly, as can be clearly seen, for instance, from 
Zi « 175 to z « 900, for the case depicted in Fig. 

Note that when the Higgs energy density starts de¬ 
creasing significantly at z > the Higgs amplitude also 
starts decreasing noticeably. However, while the Higgs 
energy density eventually stops decaying and saturates 
to an almost constant value, the amplitude \h\, instead, 
continues decreasing during a much longer time. The 
long-lasting decay of the Higgs amplitude induces the 
decrease of the potential energy of the Higgs, even long 
after the dominant energy components of the Higgs have 


saturated, as can be clearly appreciated in Fig. This 
is simply due to the fact that the Higgs keeps on os¬ 
cillating and hence transferring energy into the x fields. 
Since soon after z = Zi the energy in the Higgs becomes 
smaller than the energy in the x fields, the continuous 
transfer of energy represents only a marginal fraction of 
the energy already accumulated in the latter. Hence, the 
amplitude reached by the gradient and kinetic energy 
terms E^^Eq is not affected anymore, whereas the am¬ 
plitude of the Higgs potential energy continues decreas¬ 
ing. Eventually, the transfer of energy from the Higgs 
becomes inefficient and Ey also saturates to an approxi¬ 
mately constant value. By then, however, the Higgs po¬ 
tential energy is completely irrelevant compared to the 
gradient and kinetic counterparts. 

A very relevant aspect to note is that when all the en¬ 
ergy contributions stop growing or decreasing abruptly 
(with the exception of the Higgs potential energy, which 
keeps on falling for a long time), the energy components 
reach equipartition. In particular, some time at z > Zi, 
the kinetic energy E‘^ of the Higgs becomes equal to 
the sum of the Higgs gradient energy plus the interac¬ 
tion energy, Eq + Ei^t', see the lower panel of Fig. ^ In 
other words, equipartition in the Higgs sector holds® as 
= Eq + Aint. Similarly, in the x fields, the sum of 
their gradient energy plus the interaction term, ’equipar- 
titionates’ with their kinetic energy, E^ = Eq + Ei^t, as 
can also be well appreciated in the lower panel of Fig. 

All features described so far are, of course, not spe¬ 
cific to the particular case q = 8, (3 = 0.01 and RD, 
shown in Fig. A similar behavior is observed in the 
outcome of the field distribution for other choices of /?, 
q and w. That is, there is always initially a plateau-\ike 
stage during which the Higgs (conformal) amplitude re¬ 
mains almost constant (or changes only marginally) for 
a few oscillations. Then, the amplitude decreases fast 
when the backreaction from the x fields becomes notice¬ 
able, which causes at the same time the ceasing of the 
exponential growth of the x fields energy density. Even¬ 
tually, all fields relax into a stationary distribution with 
exact equipartition = EQ+Ei^t and E^ = EQ+Ei^t- 
On the other hand, Ey becomes completely negligible as 
compared to any other energy term E^, in correspondence 
with the decay of the Higgs amplitude, which carries on 
after equipartition is set. The duration of the different 
stages, for a given expansion rate, is directly related to 
the specific values of the parameters /? and q. In particu¬ 
lar, the duration of the initial plateau is directly depen¬ 
dent on the band structure of the Lame equation. 

We have characterized the dependence of Zi with the 
resonant parameter q; see Fig. Let us recall that Zi 
corresponds to the moment when the energy transferred 
into the x fields is sufficiently large so that the Higgs 


® In reality, it should be E'^ = E'q + Ei^t + Ey , but Ey is so small 
by then, that it does not make a difference to add it or not. 
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amplitude and energy density starts to decrease. There¬ 
fore, this is the moment that should be compared to the 


analytical estimate Eq. (37) of the Higgs decay time Zes, 


derived in Section UlI Al 

The Zi{q) behavior can be characterized by 


Zt{q) 


160 


, q G Resonant Band 


869 — 92 log q , q ^ Resonant Band 


(56) 


If a given q is within a resonant band, Zi{q) is almost in¬ 
dependent of 9 , as appreciated in Fig.[^ For RD and (3 = 
0.01, our analytical estimate Eq. ( [^ predicts z^e — 200, 
which is reasonably similar to the fit found from our nu¬ 
merical outcome, Zi{q) « 160. The analytical estimates 
are only an approximation to the real dynamics, and one 
cannot expect anything more than a reasonable order-of- 
magnitude prediction, as is indeed the case. More impor¬ 
tantly, the analytical calculation predicts that z^s should 
be only dependent on q logarithmically, which in practice 
implies that for mildly broad resonance parameters as the 
ones we have, q ~ 0(10) — O(IO^), Zes is essentially in¬ 
dependent of q, as is indeed well appreciated in Fig. 

The dependence of Zi{q) with g’s outside resonance 
bands is also logarithmic, though with a big coefficient. 
As it can be appreciated in the upper panel of Fig. 
foi' 9 ^ 10 ^ it is a factor ^2-4 larger than the analyti¬ 


cal prediction Eq. (37), but becomes of the same order 
for q ~ 10^ — 10^, modulo a factor ~l-2. Possibly, for 
q » 10^, Zi{q) will become smaller, but as said before, 
such regime is never valid in our case of study. 

In light of the results of this section, we see that the 
Higgs decay should be identified, rather than with 
with the abrupt drop of the Higgs energy density, some 
time afterwards at z > z^. After the drop, the kinetic 
contribution E'^ (which is the dominant energy com¬ 
ponent of the Higgs) enters into a stationary regime, 
equipartitioned with Eq -|- E^^i. The onset of this regime 
signals the end of the decrease of the Higgs kinetic energy. 
We therefore provide a definition of the decay time of the 
Higgs, Ze, as the moment when equipartition (within the 
Higgs sector) holds better than a given percentage. In 
practice, we operationally determine Zg as the moment 
when the equality E^ ~ Eq + Eint holds to better than 
1%. Defining like this the Higgs decay might seem arbi¬ 
trary, but when looking carefully at the evolution of the 
energy components, we see that the end of the drop of 
the Higgs kinetic energy coincides always with the 
onset of its equipartition with Eq + E-mt, for all reso¬ 
nant parameters. From then onwards, i.e. for z > Zg, 
all energy components (with the exception of the Higgs 
potential) enter into a stationary regime, evolving very 
slowly, preserving all the time the equipartition condi¬ 
tion, E'^ ~ E^ + Eint and E^ ~ E^ + Eint- 

The dependence of the decay time scale Zg versus q is 
shown in the lower panel of Fig. A fit to this relation 
is given by 


This is valid for /3 = 0.01 and for a RD (w) background. 
As we shall explain in Section |Vlj this fit can be gener¬ 
alized to other /3 values within our range of interest, and 
to other expansion rates (characterized by the equation 
of state w), as 


-(1 + 3 .^) 


Ze{q) « 50.7/3 3(1+-) q 


,0.44 


(58) 


As we can see, the behavior of the Higgs decay time 
is actually independent of whether q is within or out¬ 
side a resonance band. More remarkably, the growth of 
Ze{q) with q is actually quite contrary to the intuition 
obtained from solving the Lame equation. In the linear 
regime z < z^, when the Lame equation is valid, we ex¬ 
pect that the larger the resonance parameter, the faster 
the transfer of energy from the Higgs to its decay prod¬ 
ucts. Such trend is clearly observed (see upper panel of 
Fig. 1^, where Zi{q) either changes only logarithmically 
or decreases with g, for parameters within or outside res¬ 
onance bands, respectively. It is, however, Zg, as ex¬ 
plained, that should be interpreted as the decay time of 
the Higgs. The behavior of Zg is set by the nonlinearities 
of the problem, as opposed to z^, which is determined by 
the linear regime. This results in a completely opposite 
trend to Zi, given the growth of Zg with q. This remark¬ 
able fact, due to the nonlinear behavior of the system, 
represents one of the most relevant results of the paper. 

To conclude the section, we will briefly describe the 
dynamics of the system in the spectral domain. During 
the initial stages, the modes that are excited correspond 
to those in the band structure of the Lame equation. 
We clearly see this for z < Zi{q) in Fig. where we plot 
both the field spectra and its occupation number 

k^Uk- We also indicate with dashed lines the resonance 
bands. As the amplitude of the modes within the reso¬ 
nance bands grows, the system becomes more and more 
nonlinear. Rescattering among modes occurs, and the 
bands become wider. Due to the coupling of the modes 
through Eqs. (44) and (45), the initial parametric reso¬ 


nance of the Xk modes within the resonance bands, excite 
at the same time Higgs modes ipk>, which then rescat¬ 
ter off other modes Xk", and so on. As a consequence, 
the power spectrum of the fields grows exponentially and 
widens, with a typical width 0 < k < O(10)A:». As we 
have discussed in detail, at late times z > Zg the fields 
enter into a stationary stage, characterized by equipar¬ 
tition and a very slow evolution of the energy densities. 
This stage is indeed associated with a turbulent regime, 
typically expected to be developed due to the nonlinear 
character of a multifield interacting system [JSl 05] (see 
also isniini)- The onset of this regime translates into 
the field distributions entering in a self-similar evolution, 
with the occupation numbers verifying a scaling law of 
the type. 


n{k, t) ^no{kt ^). 


(59) 


Zg(g) = 507g' 


0.44 


(57) 


with p < 1 and g/p ^ 1 typically, and no(fc) a univer¬ 
sal function specific to each species. We have checked 
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— 335 
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— 1302 


FIG. 9: Left: Spectra for q = 14,101 of one of the scalar fields, Xi- Right: Occupation number of Xi field. The dashed vertical 
lines in the four figures indicate the position of the corresponding band of the Lame equation. Note that the units used to 
express the momentum are different from the ones used in Fig. 


that at late times z ^ Ze, the evolution of the Higgs oc¬ 
cupation number follows quite accurately Eq. (591, with 
p « 1/4 and q/p « 2.7. The late-time evolution of the 
occupation number of the x fields, however, can be fitted 
into the form of Eq. ( |59[ ) only to some extent, since any 
value between p = 1/7 and p = 1/12 does an equally 
good job (as long as p/q changes accordingly between 3 
and 4), and the high-momentum tails are always some¬ 
what offset with respect to the no{k) tails. Eventually 
the system is expected to relax into a thermal distribu¬ 
tion. The turbulent regime is, however, not very efficient 
in transferring energy from the long-wave modes to the 
high-momentum region, so an eventual total thermaliza- 
tion is indeed a long way off from the onset of the sta¬ 
tionary regime (also from our typical running times in 
the simulations). 


In the next section, we will present a similar analysis 
of the properties of the Higgs decay process, but finally 
introducing the gauge nature of the interactions. Before 
we move on, let us recall again that all our results of 
Section EYi correspond to RD and were obtained for a 
fixed value /? = 0.01. We will devote Section VI to an 
analysis of how the results change when varying the Higgs 
initial amplitude (i.e. /3) and the background expansion 
rate (i.e. w). 


V. LATTICE SIMULATIONS, PART 2: 

ABELIAN-HIGGS MODELING 

In this section, we study the properties of the Higgs 
decay modeling the system with an Abelian-Higgs frame¬ 
work. In this approach, and in contrast with the global 
scenario, we introduce for the first time a gauge structure 
in the interactions. The differences and similarities in the 
results of these two scenarios will be scrutinized. We will 
approximate the action of the electroweak sector of the 
standard model, invariant under the local SU(2)xU(l) 
symmetry group, by a local U(l) gauge theory. This is 
justified in principle because, as we will show explicitly 
in Section |V A[ the corrections due to the non-Abelian 
nature of the SM interactions, are not expected to play 
any significant role. 

Let us note that for practical reasons, we will continue 
considering a system where the Higgs is only coupled 
to a single gauge boson, with resonance parameter q = 
This way we will be able to compare directly 
the results from the gauge theory, with those from the 
previously studied global scenario. Towards the end of 
this section, however, we will consider the real case of 
the Higgs decaying simultaneously into the three gauge 
bosons ^ W~ and Z. Only in that way we will be 
really approaching realistically the dynamics of the SM. 
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Remarkably, as we will demonstrate in Section [V B[ the 
system of three (Abelian) gauge bosons can be effectively 
mapped into a system with only one gauge field, with 
effective resonance parameter q — qz + 2qw- Thanks to 
this, we will show explicitly that any analysis carried out 
with only one gauge boson, can be used directly, after 
applying an appropriate mapping, to fully understand 
the dynamics of the system with three gauge bosons. We 
will justify a posteriori this way, the ability and utility 
of modeling the system with a single gauge boson, as 
considered so far. 

The Abelian-Higgs model with one gauge boson is de¬ 
scribed by the action S = J C d'^x, with Lagrangian 

- £ = , (60) 

where the covariant derivative and field strength are de¬ 
fined as usual. 


Df, = df,- iAf, , - d^A^ . (61) 


Here, e is the Abelian coupling strength representing the 
coupling of either one of the or Z gauge fields. As 
before, in order to mimic correctly the Higgs-gauge in¬ 
teractions, we need to take = g^/4, with = g^ or 
gl, respectively for W or Z bosons. 

Since we are working with a system invariant under 
a local t/(l) transformation, we must take consequently 
the Higgs as a complex field. In terms of its components 
we shall write it as 


£ = ;^(£i+*£ 2 ) , V?* e 91 . (62) 


From Eq. (60) we derive the following equations of motion 


(p — DiDip FZK(p = 0, (63) 

9oA)iO - diF^i.i F 2e‘^a'^{t)'3m[(p*Dfj,ip] = 0 . (64) 


As in the global scenario, it is really useful to redefine 
the spacetime and field variables. On the one hand, we 
change to the same set of dimensionless spacetime coor¬ 


dinates = ( 0 °, z®) introduced in Section IV 


z = z = FI A , 


z'- = 


( 68 ) 


On the other hand, it is also convenient to define new 
Higgs and gauge field dimensionless variables as 


hj — 


a{z) ipj 

a* 


V = . (69) 


(with j = 1,2] i = 1 , 2 ,3) where (/?* = |‘/ 3 (t*)| is the initial 
modulus of the complex Higgs field at the end of inflation. 
To distinguish between different variables, we use a dot or 
a prime to denote differentiation with respect conformal 
or natural variables ( ' = d/dt, ' = d/dz), respectively. 
From now on, all spatial derivatives will also be with 
respect the new variables, unless otherwise stated. We 
also define a dimensionless covariant derivative as 


-z - iVi. 

oz. 


With these changes, Eqs. (63)-(66) can be written as 


hi - lPiz['DiT>i{hi + ih2)] ++ hl)hi = hi — , (70) 


h2 - 3m[T>iT>i{hi + ih2)] + P'^ihl + hl)h2 = ^ 2 —, (71) 
v;'+ - dAVj = j-(z), (72) 

d^Vl = jo(z) , (73) 


where the current j)j(x) is defined by 

jf,{x) = qf3‘^3m[{hi - ih2)'D^{hi + 1 / 12 )] ■ (74) 

Finally, we also define dimensionless electric and mag¬ 
netic fields as 


As we are dealing with a gauge theory, we have a gauge 
freedom in the choice of the field components. This allows 
us to set, from now on, the condition Aq = 0. In this case, 
the EOM of the gauge fields, Eq. (64), can be written in 
terms of its components as 

Aj + djdiAi — didiAj = 26"^{t)3m[(p*Djip] , (65) 
diAi = 2e‘^a^(t)3m[p*p] . (66) 


Eq. ( 66 ) is the Gauss law, which represents a constraint 


that the solution to Eqs. (63) and (65) must preserve 


at all times. When solving these equations in a three- 
dimensional lattice, one must of course check that the 


Gauss constraint Eq. ( 66 ) (or more specifically, its equiv¬ 
alent discretized version) is indeed preserved during the 
whole evolution of the system. We also define the gauge- 
invariant electric and magnetic fields as usual 

E, = Ai , Bi = ^eijkidjAk - dkAj) . (67) 


£^ = yI = §2^ - dkV,) = ||. (75) 


In this work, we have solved the system of Eqs. (70)-(73) 


in three-dimensional lattices. More specifically, we have 
solved a gauge-invariant set of analogous equations in a 
discrete spacetime. In all simulations, we have ensured 
that the lattice analogue of the Gauss conservation law 
Eq. (73) is preserved by the time evolution of the sys¬ 


tem to the machine precision. The reader can find more 
details of our lattice formulation in Appendix [Aj 

We have considered the following initial condition for 
the homogeneous modes of the fields. F rom Eq. ([ 6 ^ , we 
have by construction |/i*| = |/i(G)| = -|- = 1 at 

the end of inflation. As long as this condition is satisfied, 
we can freely distribute this initial value between the 
components hj* = hi(t*), thanks to the symmetries of 
the model. A convenient choice is 

hi* = 1 , 


h2* = 0 . 


(76) 
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As we are evolving the system of equations from the end 
of inflation, the Higgs initial velocity must obey the slow- 
roll condition (pi{t^) = With the choice 

of Eq. (761, the slow-roll condition reads 




^ 2 * ~ 0 


(77) 


We also set the homogeneous mode of the gauge bosons 
to zero, Hi, = = 0, until the onset of the oscillations 

at z — 2^osc ■ 

The system is solved in the following way. First, for 
the times 0 < z < Zosc, we only evolve the homoge¬ 
neous Higgs field with Eqs. (70) and ( [7l| ), while the ho¬ 
mogeneous gauge fields are kept to zero. At z = Zosc, 
we add fluctuations on top of the homogeneous modes 
of the different fields, allowing the gauge boson produc¬ 
tion to take place. Over the homogeneous mode of each 
Higgs component, we add Fourier modes with a spec¬ 
trum (l/fep) = 2o?uik physical variables), which mim¬ 
ics again the vacuum fluctuations of the ground state of 
a scalar field in a FRW background. Let us note that 
the initialization of the Higgs field given by Eqs. (76) is 


indeed crucial for justifying the fact that we ignore cross 
terms in the initial spectra of fluctuations. Thanks to 
the gauge rotation Eq. (76) and the slow-roll condition 
Eq. ([T^, we see that the two components of the Higgs are 
not mixed in the initial trajectory in the (/ii,/i 2 ) plane, 
and hence only the diagonal terms of the spectra of ini¬ 
tial fluctuations are needed. See [45l[52] for more details 
about this and other issues on the initialization of multi- 
field systems. 

Due to the gauge nature of the system, the initializa¬ 
tion of the gauge fields is more subtle and delicate than 
in the case of scalar fields. In this case, the fluctua¬ 
tions we add to the gauge fields must preserve the Gauss 
constraint Eq. (73) initially at every lattice point. Thus, 


given the spectrum of Higgs fluctuations, we fix the gauge 
fluctuations as given by the right-hand side of Eq. ([73 ). 


More specifically, we fix the gauge fields’ amplitude in 
momentum space as 


V) (/c, Zosc) — ^ ^2-^0 


(78) 


where in the lattice this is done with the corresponding 
lattice momenta (see |S3] for a discussion), correspond¬ 
ing to the choice of lattice finite difference operators that 
mimic continuous derivatives. The implementation of 
these initial conditions is described in more detail in Ap¬ 
pendix In particular, we discuss the importance of 
setting appropriately the Higgs fluctuations so that we 
can impose correctly Eq. (78). From z > Zqsc onwards. 


the Gauss law is then preserved to machine precision by 
the gauge-invariant evolution of the system. How this is 
checked is discussed in Appendix |A| 

We now present the main results of the lattice simula¬ 
tions carried out for the Abelian-Higgs model. Like in the 
global scenario of Section [TV] we have run simulations for 
the resonance parameters given in Table [^ ranging from 




FIG. 10: We show in blue the volume-average value of the con¬ 
formal Higgs field \h\ as a function of time for the resonance 
parameters q — 23 and q — 167, and in purple the maximum 
amplitude of the oscillations. The red line indicates the ap¬ 
proximate time at which the initial plateau finishes and the 
Higgs decay starts. The orange line indicates the value of \h\ 
at which this function stabilizes at long times, hf. 


q = 5 to q = 3000. These values correspond to A values 
between 2.5 x 10“^ and 1.5 x 10“^ for the W boson, and 
5 X 10“® and 3 x 10“^ for the Z boson. We have also 
run different simulations for (3 = 10“"^, 10“^, 10“^, 10“^ 
and 0.5. The justification of this choice of parameters 
has been explained in detail in Section |IV[ Again, all 
results presented in this section will be obtained for a 
RD background {w = 1/3) and for the (3 = 0.01 value. 
In Section |V1| we will explain how these results can be 
extrapolated to other values of w and (3. 

One of the main differences of the Abelian-Higgs model 
with respect to the global scenario is that now the Higgs 
field is described by a set of two components hi,h 2 , com¬ 
bined in a complex variable h = hi + This means 
that the quantity of interest that we must study is the av¬ 
erage value of the Higgs modulus, \h\ = \/h\ -\- Note 
that in the global case of Section |IV[ we analyzed the 
analogous quantity to \h\, corresponding to the absolute 
value of the real Higgs. This way, the results presented 
here about the decay of the Higgs amplitude can be easily 
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FIG. 11: We show the volume-average quantity \h\/a{t) oc ip 
for a RD universe. The dashed, red vertical line shows the 
time Zi, and the dashed, black line shows ~ l/a(t). 


compared with those of the global scenario. 

We have plotted in Fig. [^the volume-average of the 
Higgs modulus \h\ as a function of time, for the two 
resonance parameters q = 23 and q = 167. We have 
plotted the corresponding oscillations’ envelope curve by 
joining all local maxima with a smooth line. Remem¬ 
ber that, according to what we discussed in Section m 
all resonance parameters can be divided into two groups: 
those placed within a resonance band of the Lame equa¬ 
tion, with an interval of excited momenta of the type 
0 < fc < ^maxj and those which have a smaller band 
of the type 0 < fcmin ^ k < fcmax- We recall that ex¬ 
amples of these two groups are shown with continuous 
and dashed lines, respectively, in Fig. In this regard, 
q = 23 belongs to the second group, and q = 167 to 
the first. The initial period of oscillations, before the 
amplitude of the Higgs drops significantly, fits well the 
analytical estimate of Eq. (15). This is expected even in 
the present case with a complex field, since before the 
Higgs notices the presence of the gauge fields, the dy¬ 
namics of the Higgs radius is still effectively equivalent 
to the absolute value of a real degree of freedom. When 
the Higgs amplitude starts decreasing due to its transfer 
of energy into the gauge bosons, the period of oscillations 
is slightly modulated, but never significantly. 

We find that the Higgs amplitude behaves qualita¬ 
tively in a similar way as in the global scenario. This 
can be rapidly seen by comparing Fig. [T^ to the equiv¬ 
alent Fig. of the global scenario. In both scenarios, 
there is first a stage of few oscillations during which 
the (conformal) Higgs amplitude does not decay, corre¬ 
sponding to a plateau in the envelope function. After 
that, at times 0 > Zi{q), the Higgs amplitude starts de¬ 
caying strongly. This time is indicated in both panels 
of Fig. with a red dashed vertical line. Finally, the 
rescaled Higgs amplitude approaches a constant value at 
late times, \h\ —>■ /i/, which is indicated in both figures 
with an orange dashed horizontal line. It is important 


to emphasize again that the plateau is only manifest for 
the conformal Higgs field \h\, since the physical Higgs 
field \p\/p^ decays as ex l/a(t), due to the expansion of 
the Universe. The key observation is that, for z < Zi(g), 
it decays as the inverse of the scale factor cc 1/a, while 
for z > Zi{q) the decay is much faster due to the energy 
transfer to the gauge fields. This can be clearly seen in 
Fig.[TT] shown for q = 23. Therefore, we conclude that 
the qualitative behavior of the system is very similar, 
almost identical, to the global scenario. 

The time scale Zi{q) signals, as in the global model¬ 
ing, the moment at which the decay products (in this 
case, gauge bosons) have accumulated sufficient energy 
to start affecting the dynamics of the Higgs condensate. 
As before, this is understood better if we plot the dif¬ 
ferent contributions to the energy, as a function of time. 
The energy density of the Abelian-Higgs model Eq. (60) 
is found to be 


p[z) = 




Hz) 


Et(z), 


K = 


(79) 


where U* is the value of the Higgs potential at the end 
of inflation. The function Et{z) is formed by the sum of 
the following contributions: 


Et{z) = Ak + Agd + Ee + Em + Ey . (80) 


Here Ek and Ey are the kinetic and potential energies 
of the Higgs held 


TZV — 

£/K — 


Ey = 


W 2a? 



2 I ,,2n2 

T. ’ “ V'^1 + '^21 


Egd is a gauge-invariant term formed by the product of 
two covariant derivatives of the Higgs held (hence con¬ 
taining the spatial Higgs gradients plus the interaction 
terms) 


Egd 


i 

ih2))*'Di{hi + ih2)\ , (81) 


and Ee and Em are the electric and magnetic energy 
densities 


Ee = 


1 

K 2e2a4 




2 


Em = 


U* 2e2a4 ‘ 






2 

i ’ 


qY ^ " ■ 


(82) 

(83) 


We have plotted in the upper panel of Fig. [T^ these 
quantities as a function of time for the resonance parame¬ 
ter q = 9, which corresponds to a case within a resonance 
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band. We also show in the lower panel of Fig. the con¬ 
tribution of each energy component to the total, Ei/Et^ 
removing the oscillations of each component, and hence 
showing only the corresponding envelope functions. We 
see that initially the dominant contributions come from 
the kinetic and potential energies of the Higgs field. This 
corresponds to the oscillations of the condensate around 
the minimum of its potential, before it ‘feels’ the gauge 
fields. Meanwhile, the other components of the energy, 
E'Ei Em and E’en, grow really fast, due to the energy 
transfer from the Higgs into the gauge fields. Note that 
for the whole evolution of the system (until equipartition 
is reached), the electric energy clearly dominates over the 
magnetic energy. 

As in the global analogue, although gauge bosons are 
being strongly created, the Higgs condensate is at first 
unaffected. At z « Zi{q) (indicated by a dashed red verti¬ 
cal line in the figures) the gauge energy has grown enough 
to start affecting significantly the Higgs condensate, and 
a sharp decrease of both the Higgs potential and kinetic 
energy start from then on. Physically, this happens when 
the fraction 6 = E-^/Et < 1 becomes sizeable. In par¬ 
ticular, like in the global scenario, when d > 0 . 1 , we can 
clearly see the correspondence between the backreaction 
of the gauge fields over the Higgs field and the decrease 
in the Higgs amplitude. 

As in the global scenario, for z > Zi{q) the Higgs ki¬ 
netic and potential energies decrease sharply. The po¬ 
tential energy very soon becomes irrelevant compared to 
the other energy contributions, while the kinetic energy 
approaches an almost constant amplitude. Simultane¬ 
ously, Egb and E^, stop their growth, and also saturate 
to almost constant values. However, the magnetic energy 
continues to grow even after Eqt and Ee have been sta¬ 
bilized. Finally, at z = Ze, the system arrives again at a 
stationary regime, in which equipartition between differ¬ 
ent components is clearly achieved. In this regime, 30% 
of the total energy goes to the Higgs kinetic part, 30% to 
Egd, 20 % to electric energy Ee, and 20 % to magnetic 
energy Em- The potential energy Ey also saturates to 
a constant, but it is very subdominant with respect to 
the other contributions. Quite remarkably, these numer¬ 
ical percentages are independent of the values q and /3 
taken in our simulations. In other words, the final frac¬ 
tions of energies are universal within the Abelian-Higgs 
formulation®. 

Let us analyze the equipartition regime in the gauge 
scenario in more detail. We observe that the kinetic en¬ 
ergy of the Higgs field Ek eventually becomes equal to 
Ey + Egb- Since Egb is gauge invariant, it contains both 
the Higgs gradient terms plus the Higgs interactions with 
the gauge fields. One can then naturally identify this 


We expect this to be also the case if the non-Abelian nature of 
the interactions was considered, but only simulations of the full 
SU{2) X U(l) gauge group of the SM sector, can really prove it. 


quantity with the analogous combination in the global 
scenario, given by the sum of the interaction term plus 
the Higgs gradients, E^^ + E^. 

Similarly also to the global scenario, the potential en¬ 
ergy keeps decaying even after equipartition has been 
established. In principle, we could then think of using 
just the equipartition relation E^ — Eqe, neglecting the 
contribution from the potential energy, as we did in the 
global case. However, in the moment when the rest of 
energy contributions are stabilized, the potential energy 
still represents ~ 1% of the total. So, although the po¬ 
tential energy becomes eventually subdominant, its rate 
of decay is slower than in the global scenario. This per¬ 
centage, although small, is still significant at the moment 
in which equipartition is achieved. Therefore, it is better 
to follow the equipartition condition Ek — Ey + Egb- 
The evolution of the different energy components and the 
achievement of equipartition can be well appreciated in 
Fig.[g 


It is useful to define the Higgs decay time as the mo¬ 
ment when the Higgs kinetic energy (which dominates 
over the potential energy) results stabilized at the on¬ 
set of the stationary regime. As in the global scenario, 
we will call this quantity Ze{q)- Naturally, there is again 
some degree of arbitrariness in this definition. In the 
global scenario, we observed that a good operative crite¬ 
rion for defining Ze was based on the degree of equipar¬ 
tition achieved. In our present gauge context, we have 
observed that an appropriate criterion is to take the mo¬ 
ment when the relative difference between Ek and the 
sum Egb + Ey becomes less than 1%. We have indi¬ 
cated this time in Fig.[^ with a dashed vertical line. As 
we just mentioned, in the global scenario we did not con¬ 
sider the contribution Ey of the Higgs potential energy 
into the equipartition equalities, simply because its con¬ 
tribution was already irrelevant when equipartition was 
reached. However, in the Abelian-Higgs scenario, the ad¬ 
dition of this contribution to the covariant-gradient one 
Egb is crucial. Even though Ey is also marginal in this 
case, if we were to consider just Egb in the equipartition 
analysis, it would achieve equipartition with Ek (say to 
better than 1%) long after the Higgs kinetic energy den¬ 


sity has started to saturate. As we can observe in Fig. 12 


our criterion Ek — Egb + Ey holding better than 1%, 
coincides very well with the moment when all relevant 
energy densities have just stopped either growing or de¬ 
creasing. Hence it defines very well what we mean by the 
end of the Higgs decay. 

We have characterized again the dependence of Zi and 
Ze with the different q’s considered. We show in the upper 
panel of Fig. 13 the behavior of Zi[q). In the figure, blue 
squares correspond to q values within a resonance band, 
and purple circles correspond to values outside bands. 
We see a clear trend, such that simulations with q within 
resonance bands have a smaller Zi{q) than those with 
q between adjacent bands. Like in the global scenario, 
the order of magnitude of Zi for blue squares is approxi¬ 
mated quite well with the analytical estimate Zeff ~ 200, 
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FIG. 12: Upper panel: We plot the different contributions to 
the total energy of the system as a function of time, EijEt 
[see Eq. (80l], for g = 9. All functions are oscillating, so we 
take the envelope of the corresponding oscillations for clarity. 
The dashed vertical line signals the Higgs decay time Zi{q). 
Lower panel: We plot the same quantities with the same color 
code as in the upper panel, but now Equ and Ey appear 
dashed, and we have added a new pink line corresponding 
to Eqo + Ey, which is the quantity that equipartitionates 
with Ek- Let us note that equipartition in the gauge sector, 
between the electric and magnetic contributions, is achieved 
later than in the scalar sector, at some time 2 : > Ze{q). 


FIG. 13: Top: Different values of Zi{q) obtained for different 
resonance parameters q, for a RD universe and for (3 = 0.01. 
Blue squares correspond to q values that are within a reso¬ 
nance band of the Lame equation, while purple diamonds are 
points which are not. The purple line corresponds to the best 
fit (841, while the dashed blue line corresponds to the analyt¬ 
ical estimate ZeS ~ 200, obtained from Eq. (37l (ftj, = 0.2). 


Bottom: Red points indicate the obtained Higgs decay times 
2e(g) as a function of q, for the same Abelian-Higgs simula¬ 
tions, while the red thick line shows the best fit ( |85[ ). The 
dashed yellow line shows the best fit of this same quantity 
obtained from the global simulations in Eq. (57l. 


obtained from Eq. (37), with = 0.2. At the same time, 
the purple circles can be fitted as 


^ 1066 — 127 log g , g ^ Resonant Band, (84) 

but their dispersion is much worse than in the global case 
(recall the top panel of Fig. [^. 

In the lower panel of Fig. we also plot Ze as a func¬ 
tion of the resonance parameter g. We have obtained the 
following phenomenological fit 

Ze(g) = 588g°'42 ^ (85) 

indicated in the figure with a red continuous line. Note 
that we have plotted as well the corresponding fit ob¬ 
tained from the global simulations, Eq. ( |^ , with a 
dashed line. Both fits coincide pretty well, indicating 
that the Higgs decay time Ze{q) obtained in the global 


scenario constitutes already a very good estimation. To 
some extent this is surprising, since one could expect that 
the extra terms in the gauge field’s EOM could play some 
role, like for example modulating the decay time Ze{q) dif¬ 
ferently than in the case of only scalar fields. However, 
our results prove that this is not the case. In fact, they 
imply that the interaction term g^A^A^ifP' (which is the 
only one kept in the global scenario) is the most relevant 
one when determining the Higgs decay time scale and the 
onset of the stationary regime. 


Let us note again that the fit Eq. (85) is only valid for 


(3 = 0.01 and for a RD background. Using the theoretical 
extrapolation that we will present in Section [Vij this can 
be generalized to other /3 and w values as 


-(1 + 3 .^) 


Ze{q) « 58.8/3 3(1+-) q 


,0.42 


( 86 ) 
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FIG. 14: Electric spectra k^\Sk.\’^ and magnetic spectra k^\Bk,\^ for different times and for g = 5 (upper figure) and g = 9 (lower 
figure). The dashed, vertical lines indicate the corresponding position of the resonance band. The corresponding times at which 
the spectra are plotted are written at the right. Note that the units used to express the momentum are different to the ones 
used in Fig. 


An alternative source of information about the 
Abelian-Higgs system comes from the spectra of the dif¬ 
ferent fields. Since we are dealing with a gauge theory, 
all quantities of physical interest must be gauge invari¬ 
ant. We then plot in Fig. [^the spectra of the electric 
and magnetic fields k^\Sk\'^ and k^\Bk\‘^ at different times. 
In order to see the dependence of the spectra evolution 
on the analytical properties of the Lame equation, we 
plot both spectra for two different resonance parameters, 
g = 5 and q = 9. The latter is placed in the middle of 
a resonance band, while the former is between the first 
and second resonance bands. The dashed vertical lines 
in the figures indicate the location of the respective res¬ 
onance bands. In the case g = 5, one can clearly see 
that both spectra grow with time, as a consequence of 
the resonant excitation of gauge bosons. At initial times, 
there clearly appears a peak in both spectra, centered 
in the corresponding main resonance band. This con¬ 
firms that the behavior derived from the Lame equation 
describes well enough the real dynamics during the ini¬ 
tial stages, even for the gauge theory. When the gauge 
bosons start to affect significantly the Higgs condensate, 
i.e. for 2 > Zi{q), both spectra start to displace to the 
right, populating modes of higher momenta. In this pro¬ 
cess, new subdominant peaks appear. As time goes on. 


the peaks disappear, and when the Higgs condensate has 
decayed [i.e for z > Ze(g)], the stationary state is estab¬ 
lished. For the case g = 8, the time scale Zi(g) « 150 
is much smaller than for g = 5, and the resonance band 
is much wider. This is expected, as we include modes 
down to A: = 0. In this case, we see that the population 
of higher modes is much faster than for g = 5, and we do 
not observe additional subdominant peaks in the spectra. 

As a final remark, let us note that in the gauge sce¬ 
nario, none of the field spectra could be well fitted with 
similar scaling laws to those of the turbulent regime in the 
global case. After equipartition is reached in the gauge 
scenario, the field distributions evolve smoothly, slowing 
transferring power into higher modes, pretty much like in 
the global case. However, the evolution towards equilib¬ 
rium cannot be really grasped by simple fitting formulas 
like Eq. pi). 


A. Beyond the Abelian-Higgs 


The real nature of the SM interactions is non-Abelian, 
since the EW sector of the SM is SU (2) x U (1) gauge in¬ 
variant. In the EOM of the gauge bosons there are there- 
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fore nonlinear terms^ of the form ~ g^A^,gAdA,gdA^, 
where we omit charge and Lorentz indices for simplicity. 
Following m, one obtains that within the Hartree ap¬ 
proximation, the terms ~ gAdA, gdA^ vanish, so that in 
principle only the terms ~ g^A^ contribute effectively to 
the dynamics of the gauge fields. We can write the effec¬ 
tive mass entering into the gauge fields’ EOM, as given 
by their interactions with the Higgs, plus a contribution 
from their own non-Abelian self-interactions. Symboli¬ 
cally, we will write this as 

m-A = • (87) 

The Abelian-Higgs simulations capture the first term 
g'^ip'^, which is due to the interaction with the Higgs, and 
is responsible for the resonant excitation of the gauge 
fields. The self-induced mass due to the gauge-field self¬ 
interactions is, of course, not present in the Abelian ap¬ 
proach. This second term describes the nonlinearities of 
the non-Abelian nature of the SM interactions. Hence, 
only when the gauge fields have been excited with a suf¬ 
ficiently high amplitude (A^) > g'^ip'^ may their presence 
have any relevance. The question, then, is when do the 
gauge fields reach the critical amplitude A ^ Ac = gifl 
The answer can be easily found by analyzing the effec¬ 
tive mass of the Higgs. The non-Abelian nature of the 
interactions does not add any extra contribution into the 
effective mass of the Higgs field, given by 

ml = X^^+g^{A^) . (88) 

These terms are already captured in our simulations, so 
the only difference in a non-Abelian simulation would 
come from the fact that A^ is affected by the nonlin¬ 
earities of its own EOM. The gauge fields backreact into 
the Higgs dynamics at the time z = Zi{q), which cor¬ 
responds physically with the moment when the ampli¬ 
tude of the gauge fields has grown - due to parametric 
resonance - up to (A^) > Xq?jg^. This condition cor¬ 
responds, however, to a typical amplitude of the gauge 
fields A ^ A(zi) = -sfXLpIg, which is much smaller than 
Ac. In particular, for the typical broad 

resonant parameters q ~ 0(10) — O(IO^). The effective 
mass of the gauge bosons at z « is 

mA(zi) = « 5^^ (^1 + , (89) 

where ^ ^ 1- It is then clear that m\{zi) « as if 

there were no effect from the gauge-field self-interactions. 
By the time the gauge-field resonant production backre- 
acts on the Higgs dynamics, the gauge fields stop grow¬ 
ing, as explained in detail in Section |V] Therefore, the 


^ For the sake of clarity of the physics, we switch back to physical 
variables in the discussion of this subsection. 


non-Abelian terms (neglected in the Abelian-Higgs ap¬ 
proach), are not expected to play any significant role 
in the dynamics of the system.® It is, however, likely 
that the presence of the non-Abelian terms will possibly 
change the details of the achievement of the equiparti- 
tion regime. Therefore, although we do not expect the 
time scale Zi{q) to change, the time scale Zf,(q) might per¬ 
haps change moderately in the presence of non-Abelian 
corrections. However, only non-Abelian lattice simula¬ 
tions, beyond our present work, can really quantify these 
questions. 

In light of this analysis, we see a posteriori that ne¬ 
glecting the nonlinearities due to the non-Abelian nature 
of the SM interactions was well justified. 


B. Abelian-Higgs model with three gauge fields 

So far, we have studied the postinflationary Higgs dy¬ 
namics in the lattice, mimicking its interaction with a sin¬ 
gle gauge boson using an Abelian-Higgs modeling. This 
has allowed us to obtain a bunch of interesting results, 
which depend greatly on the choice of the gauge boson 
resonance parameter, q = (/^/(4A), with g"^ being the 
corresponding standard model coupling of either W oi Z 
bosons. Naturally, we should include the three massive 
gauge bosons in our simulations (i.e. the IE®", W~ and 
Z), as in the EW sector of the standard model. Remark¬ 
ably, the results presented so far for a single gauge field 
can be easily translated into the three-boson case, with 
an appropriate field redefinition. We explain this in what 
follows. 

In the case of a Higgs decaying into three Abelian 
gauge fields, the Higgs equation can be written as 

a" 

h"+P'^\h\'^h = h— (90) 

a 

where h = hi -\- ih 2 ^ and the covariant derivative is now 

+ W,- + Z,) . (91) 

Here, W+, IE", and are the corresponding fields of 
the IE®", IE“, and Z bosons, respectively. We describe 
the three fields in the temporal gauge, so that their 0 
components are null. The EOMs of either of the IE 
bosons are then 

IE" + djd^W^ - d,d^Wj = qwl3‘^Am[h*V,h] , (92) 
diW' = qwP‘^'3m[h*h'] , (93) 


It is possible though that for the mildest broad resonance pa¬ 
rameters, such as g ~ 0(10), there might be some effect from 
the non-Abelian terms, since in this case the correction in 

j—I ^ ^ 

Eq. (|89|l is only marginally smaller than unity. 
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with qw = 5^/(4A). Equivalently, the EOMs of the Z 
boson are 


Z” + djdiZi - dAZj = qzl3‘^'3m[h*'Dih] , (94) 

d,Zl = qzP^'3m[h*h'] , (95) 

with qz = 5|/(4A). Note that there is a Gauss law for 
each gauge field, representing as before, dynamical con¬ 
straints of the system. Interestingly, this system can be 
reduced, with an appropriate redefinition of the gauge 
fields, to the case of a Higgs decaying into a single gauge 
field studied above. To see this, let us define the following 
effective gauge field 


Sf, = W+ + W- + Z^, 

and the resonance parameter 


q = qz + 2gw = 


g| + ‘^9w 

4A 


If we consider the mapping 




(96) 


(97) 


(98) 


automatically S'o = 0, and we can then reduce both the 


W EOM (92)-(931 and the Z EOM (94)-(951 to just 


S'/ + djdiSi - dASj = ql3'^3m[h*Vih] , (99) 


diS[ = ql3'^'3m[h*h'\ 


( 100 ) 


where the covariant derivative of Eq. (|91[) is now simply 


'Df, = d^- iSfj, 


( 101 ) 


Therefore, the three gauge bosons can be described® 
by a single effective gauge boson Si, coupled to the Higgs 
with the resonance parameter q of Eq. (97). This prop¬ 


erty is very useful, since we just need to introduce only 
one effective gauge field, Eq. (96), and the system is 
then fully described by Eqs. (90), (99), (100), and the 


covariant derivative (101). As an example, if we have 


qw = 14 and qz — ‘2.qw = 28, all three gauge bosons 
can be described by the EOM of a single gauge field with 
resonance parameter g = 28 -I- 14 -|- 14 = 56. In other 
words, the system behaves in such a way that the three 
gauge bosons have the same effective resonance parame¬ 
ter. From Eq. (98), we find the following relation between 


the W and Z amplitudes 


Z,(z) 


qz 

qw 


W+{z) 


qz 

qw 


W-{z) , 


( 102 ) 


^ Note that this property can be generalized to a Higgs coupled 
to N Abelian gauge bosons, Vf (c = 1, 2, • • • N), with different 
resonance parameters Qq. If we define Si so that Vf = (qc/q)Si, 
all fields can be described with the same eom of a single field 
with effective resonance parameter q = qc. 


which at very high energies, when qz ~ 2 qw , reduces 
simply to Zi{z) « 2W/'(z) « 2W~(z). Eq. (102) follows 


in all spacetime (and in the lattice, in all sites at all 
times). 

We have just seen that the dynamical equations of the 
Higgs coupled to three gauge bosons can be reduced to a 
system with the Higgs coupled to only one gauge boson, 
with resonance parameter q = qz + 2qw- The equiva¬ 
lence between these two systems is actually a mathemat¬ 
ical identity. For the sake of verification, we have checked 
that the results in terms of Zi and Zg, are indeed identi¬ 
cal when comparing the simulations of one effective gauge 
boson S/^ with the resonance parameter q = qz + 2qw, 
and the simulations of two W bosons with the resonance 
parameter qw each, plus a Z boson with resonance pa¬ 
rameter qz- 

Given the above equivalence, in principle, we should 
then be able to translate the analysis from the simula¬ 
tions discussed so far, for only one gauge field, into the 
real scenario with the three gauge fields W^, Z. Strictly 
speaking, however, both scenarios are not really identical, 
if we compare them for the same q and /3. To understand 
this, let us identify the gauge boson of the single gauge 
field simulations presented so far, e.g. with one of the 
W bosons. For a fixed value of its resonance parame¬ 
ter q = g'^/AX, we conclude that the Higgs self-coupling 
is A = XiB = q’wl'^q- In the case of the three gauge 
bosons, however, the effective field 5^ (exactly equiva¬ 
lent to the three-gauge-boson system), has a resonance 
parameter q = {f2g^ -l-^D/dA. From there we deduce 
that for the same q, the Higgs self-coupling in this case 
is A = Ass = (2^^ -I- 5 |)/4(7, which differs in a factor 
(2 -I- {gz/gwY) with respect to Ais- In other words, we 
would be comparing systems with different Higgs self¬ 
couplings, and hence not equivalent. Since we want to 
compare the two systems for the same j3 = VXa, the dif¬ 
ference in A translates into a different a, and hence into 
different initial Higgs modes. According to the initial 
spectra given by Eqs. (B7), (B8), the fluctuations de¬ 


pend explicitly on A, as a reflection of their dependence 
on a (after having fixed /3). Given our choice of variables, 
the fluctuations added at the time Zosc, are then slightly 
different in the two scenarios. As a consequence, there 
are also differences in the gauge initial fluctuations, as 
we impose, following the procedure of Eq. (78), 


S'i{k,Zosc) = i^joik,Zosc) , (103) 


with jo{k,Zosc) the Fourier transform of jg^ZjZ) = 
qesl3^'3m[h*h'], evaluated at z = Zosc- 

It is crucial, then, that we figure out the importance 
of these differences. If they are irrelevant, we can then 
simply use the results presented so far for a single gauge 
field, for describing the real case of three gauge bosons. 
In order to find this out, we have compared the results for 
Zi and Ze from simulations with only one gauge boson, 
identical q and /3, but different A ( = Ait, and Aat,, ac¬ 
cording to the discussion above). In the new simulations 
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with one effective gauge boson S^, we have observed the 
same dynamics as in the case of only one W boson. We 
first observe a stage in which the volume-averaged Higgs 
amplitude \h\ shows a. plateau. In this regime, as before, 
the gauge energies grow very fast, but their contribution 
is still not important enough to affect the Higgs conden¬ 
sate. The times Zi at which the plateau ends are reduced 
slightly with respect to the W-boson case when q is out¬ 
side a resonance band, but they are almost identical when 
it is within a band; see Fig. |15[ There are, however, vir¬ 
tually no differences in the time scale Ze, which signals 
again both the end of the Higgs decay and the onset of 
equipartition. The new fit of Ze from the simulations with 
an effective gauge boson is 


Ze{q) = 58lq^-^^ = 58liqz + 2qw) 


0.42 


(104) 


very similar to the old fit Eq. (85). Note that, as men¬ 


tioned before, this fit is done for a RD universe with 
/3 = 0.01. Anticipating again the results that we will ex¬ 
plain in Section [VTl the generalization of this fit to other 
13 values and expansion rates (characterized by w) is 


zM ~ 58.1/3^^ {2qw + qzf-^^ , (105) 


This equation probably represents the most relevant re¬ 
sult of our paper. We see that the real decay time Ze 
of the Higgs into the three gauge bosons W'^,Z is, us¬ 
ing again the approximate high-energy relation qz « 
2qw, a factor {{qz + 2qw)/qw)^'"^^ ~ ~ 1-79 

times longer than if we only considered the decay of 
the Higgs into a single W boson [equivalently, a factor 
{{qz + 2qw)/qz)^''^^ ~ 2° "^^ « 1.34 longer if we consid¬ 
ered the decay of the Higgs into a Z boson]. 

It is perhaps worth noticing that it seems surprising at 
first glance, that the decay takes longer when the reso¬ 
nance parameter is effectively larger, q = 2qw+qz > qw', 
naively one would expect a faster decay if there are more 
bosons into which to decay. This is, however, a reflection 
again of the nonlinear behavior of the system at z ^ Zi, 
responsible for the previously discussed counterintuitive 
growth of Zf.{q) with q. 

We also expect the energy equipartition not to change 
with respect to the single W boson case, simply because 
the way in which the gauge fields are initially excited, 
should not affect the late-time dynamics when the non- 
linearities are important. As confirmed by the lattice 
simulations, this is indeed the case. We have checked 
that the hnal equipartition state is identical to the pre¬ 
viously studied case of one single boson, reaching at late 
times. 


^ « 0.3 , ^ 

Et Et 


0.3 


Ee 

Et 


0.2 


Em 


0.2 , (106) 


and Ey/Ei ^ 1. Note that, in the case of three gauge 
bosons, we have three different electric and magnetic 
fields. From the relation Zi{z) = 2Wi{z) (valid at high 
energies), and given the definition of the electric and 
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FIG. 15: Top panel: Filled points show the Zi{q) times ob¬ 
tained from simulations with an effective gauge boson 
whereas empty points show the analogous results from simu¬ 
lations with a single HA boson shown in Fig. |13[ Blue squares 
and purple diamonds correspond to q values inside and out¬ 
side a resonance band of the Lame equation. Bottom: points 
represent the Ze(g) values obtained for the effective Sfj, boson, 
whereas the blue line corresponds to the phenomenological fit 
of Eq. |T04|. 


magnetic energies (82)-(83), we see that 50% of the to¬ 
tal electric energy corresponds to the Z boson, while the 
other 50% is divided equally between the other two W 
bosons. The same distribution takes place for the mag¬ 
netic energy. 


As a final remark, let us note that before nonlinear ef¬ 
fects become important, the behavior of the three gauge 
fields is described by the same Lame equation with res¬ 
onance parameter q = 2qw + qz- Due to this, we have 
observed that for z < Zi(q), the spectra of the three fields 
excite the same range of momenta, corresponding to the 
resonance band of the Lame equation for such resonance 
parameter. This is important, as one could naively think 
that the spectra of and Z are independent, with 
different rangesof excited momenta accordingly to their 
different resonance parameters qw and qz- On the con¬ 
trary, the introduction of three gauge fields in the system 
makes them evolve, as seen, as a single effective gauge bo- 
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son, with the same effective resonance parameter given 
by Eq. 


VI. VARYING THE HIGGS INITIAL 
AMPLITUDE AND THE EXPANSION RATE 

All results from sections|lV]and|V]have been presented 
for a scale factor evolving in a RD universe (w = 1/3), 
and for fi = 0.01. Naturally, in order to fully understand 
the dynamical properties of the Higgs decay after infla¬ 
tion, we have explored other /3 parameters, and we have 
also considered other expansion rates such as MD {oj = 0) 
or KD (w = 1). Fortunately, one can easily extrapolate 
the results from one particular set of parameters, say (/3i, 
wi), to another set (/32, W2), using the analytical prop¬ 
erties of the Higgs equation described in Section |Hj In 
other words, from the results obtained for (/3i, wi), we 
can obtain a very good approximation to the ones for 
(^2, W2). 


More specifically, we saw in Eq. (181 that in the case 


of no coupling to the gauge bosons, the conformal pe¬ 
riod Zt and the value of the transformed Higgs field 
at the first maximum h{zM), can be approximated as 

Zt = ci/3 3(1+0.) and /i{zm) = £ 2/? 3(1+0,) ^ -where ci and 
C2 are constants independent of lo and /3. From these 
properties we can see that, if for a given set of values 
(a;i,/3i), the volume-averaged Higgs field takes the value 
h{l3i,oji) at the time z(/3i, wi), then for (0)2,/32) the Higgs 
field at the time 

-(1 + 30,2) (1+30,1 ) 


Z{P2,UJ2) ^ 

ziPuLO,) , 

(107) 

should take the value 

-2 

2 


h{P2,i02)-p!™ 


(108) 


Notably, this property is maintained quite well even in 
the presence of a Higgs coupling to its decay products 
(either scalars in the global simulations or gauge bosons 
in the Abelian-Higgs simulations). This extrapolation is 
therefore very powerful. In Fig. we have plotted the 
volume-averaged value of \h\ as a function of time, for 
both global (top figures) and Abelian-Higgs simulations 
(bottom figures). Let us focus for instance on the top-left 
figure. We have obtained for q = 8 the behavior of \h\ 
as a function of time for j3 = 10“'*, 10“^, 10“^, 10“*, and 
0.5, directly from the simulations. Using the outcome 
from these simulations with different f3 parameters, we 
have then inverted Eqs. (1071 and (1081, and obtained 


the (extrapolated) behavior corresponding to j3 = 0.01. 
These are different predictions for the Higgs decay when 
(3 = 0.01, but obtained from the real data from simula¬ 
tions with different j3 values. We see that the four dif¬ 
ferent extrapolated theoretical predictions obtained for 
P = 10“^, 10“^, 10“* and 0.5 coincide very well with the 
real simulation for /3 = 0.01. 

The same is done in the top-right figure, but changing 
the scale factor instead of P (which we fix in this figure 


ds P = 0.01). There, we compare the result of the Higgs 
decay for a RD universe, on one hand obtained directly 
from simulations with lo = 1/3, and on the other hand 
from the corresponding extrapolated predictions from the 
lattice simulations with (u = 0 (MD) and oj = 1 (KD). 
The three lines also coincide very well. The same anal¬ 
ysis is repeated for Abelian-Higgs simulations in the two 
bottom figures, with identical conclusions. 

This property allows us to extrapolate easily the results 
for the Higgs decay time for a RD universe with P = 
0.01, presented in the last two sections, to another set 
of {ui,P) parameters. In particular, from Eq. (57) we 
obtain Eq . ((^ , from Eq. (85 ) we obtain Eq. (86), and 
from Eq. (1 1041) we obtain Eq. (105). 


VH. SUMMARY AND DISCUSSION 

The recent measurements of the Higgs boson mass [U 
[5] imply a relatively slow rise of its effective potential 
at high energies. In the regime where the EW vacuum 
is stable with the Higgs self-coupling kept positive, the 
Higgs develops a large VEV during inflation, representing 
a classical condensate, homogeneous over scales exponen¬ 
tially larger than the inflationary radius l/H*. In this 
paper we have studied the relaxation of the Higgs, i.e. its 
decay, during the stages following immediately after in¬ 
flation. In reality, the origin of the VEV during inflation, 
which sets up the initial condition for the decaying pro¬ 
cess, is not particularly relevant for our study. If another 
mechanism (different than quantum fluctuations) is re¬ 
sponsible for the development of the Higgs VEV during 
inflation, our calculations and results would be equally 
applicable. The case considered in the paper, with the 
initial amplitude of the Higgs condensate dictated by the 
equilibrium distribution Eq. due to the stretching 
of its quantum vacuum fluctuations, simply serves as a 
starting and practical point, to assess the typical Higgs 
amplitudes at the end of inflation. 

The decay of the Higgs condensate during the early 
postinflationary stages constitutes an important event in 
the evolution of the Universe, which might have inter¬ 
esting cosmological consequences. In this article we have 
focused on the details of the Higgs decay process itself. 
We have used different methods of progressive complex¬ 
ity, accuracy and proximity to the real case of the SM. We 
have modeled the SM interactions in a two-step manner. 
First, considering a global scenario, ignoring the gauge 
structure of the SM, representing the gauge fields as a 
collection of scalar fields appropriately coupled to the 
Higgs. Secondly, we have considered an Abelian gauge 
scenario, with the gauge fields and the Higgs embedded 
within an Abelian-Higgs framework, ignoring the nonlin¬ 
earities due to the truly non-Abelian nature of the SM. 
For the global model we have presented both analyti¬ 
cal (Section HI AI and lattice calculations (Section IV), 
whereas in the most precise and involved gauge mod¬ 
eling, we have just presented the outcome from lattice 


























27 



FIG. 16: We plot the volume-averaged value of the Higgs conformal field \h\ as a function of time, obtained directly from our 
simulations, for either different j3 parameters or expansion rates. Lines with the symbol ‘(r)’ have been extrapolated, using an 
inversion of Eqs. ( |107[ ) and ( |108| l, to obtain a theoretical prediction of the results of a RD universe {uj = 1/3) with /3 = 0.01. The 
top hgures correspond to global simulations with q = 8, and the bottom hgures correspond to Abelian-Higgs simulations with 
q = 6. In the left hgures, we vary /3, while in the right hgures, we vary lu. We see that the lattice results for {to, /?) = (1/3,0.01) 
coincide quite well with the different theoretical extrapolations obtained from the lattice results for other (tu, /?) parameters. 


simulations (Section [V|). 

The analytical results of the global modeling estimate 
correctly the right order of magnitude of the Higgs de¬ 
cay time. When studying such a scenario in the lat¬ 
tice, including all nonlinearities within such a scheme, 
we find that the actual Higgs decay takes longer, typi¬ 
cally a factor Ze/^eff ~ 3.17q° ‘*^ larger: see Eq. (58l for 
Ze and Eq. (37) for ZeS- This is because the analytical 


calculations are only capable of estimating the order of 
magnitude of the time scale when sufficient energy has 
been transferred into the extra scalar fields (mimicking 
the EW gauge bosons). However, that time only signals 
the moment z = Zi{q) when the Higgs condensate re¬ 
ally starts noticing that it is coupled to extra species. 
From then on, at times z > Zi{q), the Higgs energy 
density begins to decrease in a noticeable manner, be¬ 
ing transferred to the most strongly coupled species, the 
EW gauge bosons. It is this decrease of the energy of 
the Higgs that should be interpreted as the decay of the 
Higgs. Eventually, the Higgs energy density saturates 
to an approximately constant value, at some moment 
Ze{q) > Zi{q). Around the same time, the energy of the 
species coupled to the Higgs has also stopped growing. 


and saturates into slowly evolving magnitudes. 

Very interestingly, the same pattern and time scales 
are observed in the gauge scenario, though the final frac¬ 
tions of energies are different. The time scale Ze{q) that 
characterizes the end of the Higgs decay in the gauge 
case is given by Eq. (86), which represents a factor 
ZelzeG ~ 3.68q°-‘^^ larger than the analytical prediction 
Zefi of Eq. (37). We see therefore that, at the end, the 


differences between the global and gauge modelings are 
not so relevant, at least in terms of the estimation of the 
Higgs decay time Ze{q)- It is worth stressing that Ze{q) 
grows with q (both in the global and gauge scenarios), 
which could be thought as being a counter-intuitive fact. 
This is due to the nonlinearities characteristic of the sys¬ 
tem, which become relevant from z > Zi onwards. 

One of our more interesting results is the extrapolation 
laws Eqs. (107),(108). We have seen that the dynamics 
of the system depend basically on three parameters: q, 
/3, an d th e expa nding background equation of state ui. 
Eqs. (107),(108) allow us to extrapolate the lattice re¬ 


sults for parameters (wi, /3i) into a very good approxima¬ 
tion to the results of another set of parameters {u} 2 ,132)■ 
This technique works very well indeed for both global 
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and Abelian-Higgs simulations (see Fig. 16). This hap¬ 


pens because the properties of Eqs. (16) and (18) derived 


in Section III, also hold quite well in the presence of a 
coupling of the Higgs to its decay products. This has led 
us to obtain the generic formula for the Higgs decay time 
0e, Eq. (105), as a function of (3, q and w. 


Remarkably, we have also shown that the case of the 
SM, where the Higgs is coupled simultaneously to the 
three EW gauge bosons IF+, W~ and Z, behaves iden¬ 
tically to the case in which the Higgs is only coupled 
to one effective gauge boson, with resonance parame¬ 
ter q = qz + 2qw- We have found that when the three 
gaug e bos ons are considered, Ze{q) = 581((7z -I- 2qw)^''^^ 
[Eq. (104)]. The decay of the Higgs takes then a factor 
(2 -I- qzlqw)^'^"^ larger than if the Higgs were coupled to 
only one W boson, or equivalently {ff + 2qw jqz'ff '^’^ times 
larger than if it were coupled to only Z gauge bosons. 
Again, this counterintuitive result is due to the nonlin¬ 
earities that dominate the system at z > 


Interestingly, at the time z « Ze{q), in both in the 
global and gauge scenarios, we see that the distributions 
of fields reach equipartition. In the global model we find 
that the kinetic energy of the Higgs becomes equal to 
the sum of the gradient energy of the Higgs plus the 
interaction with the Xi fields, ~ Eq + Ei„f This 
equality holds to better than 1 % from z > Ze onwards. 
In the gauge scenario, we find that the kinetic energy 
of the Higgs becomes equal to the sum of the covariant 
gradient energy (which includes the Higgs-gauge bosons 
interactions) plus the Higgs potential, Ek — Eqd + Ey. 
This equality also holds to better than 1% from z > Ze{q) 
onwards. At some later time z > Ze, the electric and mag¬ 
netic energy densities also reach equipartition to better 
than 1%, Ee — Em- The distribution of energy in the 
gauge scenario is actually universal, since the system al¬ 
ways reaches equipartition, with Ek — Eqd representing 
30% of the total energy, and Ae « Em representing 20% 
each. In both global and gauge scenarios, once in the 
stationary equipartitioned regime, the potential energy 
becomes gradually more and more irrelevant. 


Before we conclude, let us note that the postinfla¬ 
tionary decay of the Higgs analyzed here is very simi¬ 
lar to the analogous decay during reheating after Higgs- 
inflation [SlilHSI]. The contexts are, however, very dif¬ 
ferent. In Higgs-inflation the Higgs plays the role of the 
inflaton and dominates the energy budget of the Uni¬ 
verse, so the decay of the Higgs after inflation truly rep¬ 
resents the actual reheating of the Universe. In the case 
we have studied in this paper, the Higgs is simply a spec¬ 
tator field during inflation, and its energy density is only 
a marginal fraction of the inflationary one; see Eq. (20). 
In Higgs-inflation, a nonminimal coupling to grav¬ 
ity is required, with ^ ~ The resonance 

in both Higgs-inflation and Higgs-spectator scenarios is 

dominated by the decay into the gauge bosons IU*,Z. 

2 

The resonance parameter, however, scales as g ^ 

2 

in Higgs-inflation, versus g ~ ^ in our Higgs spectator 


scenario. Therefore, the resonance is ^ 10^\/Aooi times 
broader in Higgs-inflation than in the Higgs-spectator 
case. However, in Higgs-inflation, the nonperturbatively 
produced gauge bosons (at each Higgs zero crossing), de¬ 
cay very fast into the SM fermions via perturbative de¬ 
cays. So for around ^ 100 oscillations of the Higgs, the 
resonance is blocked in Higgs-inflation, simply because 
the occupation numbers of the gauge bosons do not pile 
up [SS] . This phenomenon is called combined preheating, 
and it is absent (or in general it is expected to be only a 
marginal effect) in the Higgs spectator scenario studied 
here, as shown in |22) . 

To conclude, let us note that our paper is intended 
to be the first one of a series, where we plan to analyze 
further the details of the Higgs decay (i) and its cosmo¬ 
logical consequences (ii). In particular, 

(i) The results obtained here have gone far beyond the 
analytical ones available in the literature [HI HO]- We 
have presented different approaches to the nonperturba- 
tive and nonlinear dynamics of the decay process. Our 
most precise results are the outcome from our simula¬ 
tions in Section |Vj corresponding to an Abelian-gauge 
model mimicking the structure of the SM interactions. 
Even though there is a good motivation to neglect the 
truly non-Abelian nature of the interactions, only lat¬ 
tice simulations which fully incorporate the non-Abelian 
SU{2) X 17(1) structure of the SM will really tell us about 
the (un)importance of the corrections due to the nonlin¬ 
earities in the gauge sector. Besides, the details of the 
stationary stage might very well (and indeed most likely 
will) change when the full non-Abelian structure of the 
SM is restored. Therefore, even if the time scales of the 
start of the Higgs decay and onset of stationary regime 
may (expectedly) not change much, the fine details can 
only be quantified in light of such non-Abelian simula¬ 
tions, which are beyond our present work. Moreover, in 
order to assess with even a higher degree of realism the 
final outcome of the energy distribution among fields, 
thermal corrections [581160] and fermions [6T1 - I63] should 
be effectively incorporated into such simulations. 

(ii) The postinflationary decay of the SM Higgs may 

have several observable consequences. The possibility has 
been recently proposed [641165] of realizing baryogenesis 
via leptogenesis, thanks to the Higgs oscillatory behav¬ 
ior. The time dependence of the Higgs condensate oscil¬ 
lations can create an effective chemical potential for the 
lepton number, which could lead to the generation of a 
lepton asymmetry in the presence of right-handed Ma- 
jorana fermions with sufficiently large masses. The elec- 
troweak sphalerons would then redistribute such asym¬ 
metry among leptons and baryons. Second, the fields 
excited from the decay of the Higgs may act as a source 
of gravitational waves The case of the charged 

fermions of the SM was considered [43] , but it is expected 
that the background of gravitational waves from the EW 
gauge bosons contributes to a much larger signal [43] . 
Besides, the fact that the Higgs is a condensate vary- 
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ing at superhorizon scales may give rise to interesting 
anisotropic effects [731 El in the amplitude of such a 
background of gravitational waves. Thirdly, it is indeed 
possible that the gauge field production that we have de¬ 
scribed in this paper could provide the necessary condi¬ 
tions for primordial magnetogenesis. Although it might 
be challenging to obtain a sufficiently large correlation 
length, it is conceivable that an inverse cascade process 
provides the appropriate mechanism for the growth of an 
initially small correlation length [751176]. Finally, if dark 
matter is a gauge singlet field coupled to the Higgs, it is 
also possible that the Higgs oscillations could produce the 
right amount of dark matter, such that its distribution 
could account for the correct relic abundance m 


Note added. - After completion of this work, the 
preprint [78| by Enqvist et al was uploaded to the ArXiv, 
presenting lattice simulations of the same process an¬ 
alyzed in this paper, but considering the non-Abelian 
structure of the standard model. Only low-resonance pa¬ 
rameters with q < 20 were considered, and for a fixed 
initial amplitude of the Higgs and postinflationary ex¬ 
pansion rate. The expected broadening of the gauge field 
spectra, due to the nonlinearities introduced by the non- 
Abelian terms, is indeed clearly observed after some time, 
as compared to the Abelian simulations. However, for the 
lowest case of g « 6, where the effects of such nonlineari¬ 
ties are expected to be maximum, only a factor ^ 2 of dif¬ 
ference in the estimation of Zi is observed, as compared to 
the analogous Abelian simulation. For higher-resonance 
parameters, the Abelian approximation becomes better 
and better, as the correction due to non-Abelian terms 
become more and more irrelevant, see Eq. (89). Besides, 
in the context of a large inflationary energy scale (close 
to its upper bound), it is rather expected that q ^ 10, 
as <7 ^ 0(10) requires an excessively large Higgs self¬ 
coupling. Therefore, we are positive that the work we 
have developed here is a very good approximation to the 
real, non-Abelian dynamics. We plan to study this issue 
in a future publication. 


Acknowledgments 


We thank Mustafa Amin for illuminating discussions 
on the initial conditions. We would like to thank Fe¬ 
dor Bezrukov for making publicly available the package 
to compute the running of the Higgs self-coupling in the 
website http;//www.inr.ac.ru/~fedor/SM/, F.T. ac¬ 
knowledges the CERN Theory Division for kind hospi¬ 
tality. This work is supported by the Research Project 
of the Spanish MINECO FPA2012-39684-C03-02 and the 
Centro de Excelencia Severo Ochoa Program SEV-2012- 
0249. P.T. is supported by the PPI-Severo Ochoa Ph.D. 
fellowship SVP-2013-067697. We acknowledge the use of 
the IFT Hydra cluster for the development of this work. 


Appendix A: Lattice formulation 

In this appendix, we provide a more detailed discussion 
of the lattice formulations for both the global and the 
Abelian-Higgs simulations. Let us start by writing the 
action for both scenarios in the continuum. In the global 
case, the continuous action Eq. (39) can be written in 


our natural variables as (from now on a* = 1) 


S = 


d^ 


^ I - (h' - "H/i)^ -I- dihdih 




rN 

4A 


e2/32 




(Al) 


Varying this action, we find the continuum EOMs of the 
system 


h"+ + e^h^X] = —h, (A2) 


A" - + qP^h'^Xj = —Xj , (A3) 


which are Eqs. (44) and (^45^ of the main text. 


We now want to write the equivalent of these equations 
in the lattice. We will work in a lattice cube of length 
L with points. We take the time step to be do and 
the lattice spacing to he di = d = L/N {i = 1, 2,3). We 
write 


Aq A+/i-A-A+/i + /32/i3 + e"h^A2 = —h, (A4) 


A-A+A, - A-A+A, + qP^h^X, = -A, , (A5) 


3 

2,2 


1 ? 
a 

where we have defined the discrete derivatives A+i)) = 
+ p) - p) and A-p = ^{p - p{n - p)). Nor¬ 
mally, one obtains the operators A~A+ from discretizing 
the continuum action ( |A1| ) and then minimizing it with 
respect the lattice field variables (which live in the lat¬ 
tice sites n). However, since we are not treating at the 
same level in the lattice the scale factor a(t) (which we 
do not discretize) and the field variables h, Xj (which are 
discretized), we prefer to proceed by simply substituting 
the continuum operators for the lattice equivalent 

Ap A+. Actually, if we indeed proceeded by discretizing 
the action and then finding the lattice EOM, we would 
of course obtain the lattice operators A“A+ on the left- 


hand sides of Eqs. (A4) and (A5). However, on the right- 

r/ 1—1 

hand side, the term —h would be more involved in dis- 
Crete derivatives of time. Since we know that such a term 
decays very fast, we have simply introduced the term 
on the right-hand side of the EOM by using a continuous 
function a{t) evaluated at the appropriate discrete times. 
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With respect to the Abelian-Higgs simulations, let us 
first write its continuous action Eq. (601 in terms of nat¬ 
ural variables (where h = hi + ih 2 ): 


d^- 


S= ; 


/32 


2 1 


-\h'-mf +\D,h\'^ 


i E ^ E ^o^ + ^ih*hy 


Aq 




(A6) 


where = d^Vi, — di,V^. Varying this action, we obtain 

the EOM in the continuum 


end of the running time the Gauss law is in fact only 
marginally broken, with Aq ^ 10“^^ — 10“^^. 

All results presented in this work have been obtained 
for N = 128 points for both global and Abelian-Higgs 
simulations. Apart from iV, we also need to fit the range 
of momenta that we want to cover in our simulations. 
This is a crucial step, as this range must be chosen care¬ 
fully in order to capture all the relevant phenomenology 
of the Higgs decay. Let us call Pmin the minimum momen¬ 
tum covered by the lattice. We can then fix the length 
of the cube L and the maximum momentum covered by 
the lattice, Pmax, in terms of N and Pmin as 


h" - D,Dih +(3^\h\^h = h—, (A7) 

V” + djdiVi - didiVj = ji{z) , (A8) 

= Joiz) , (A9) 

where the current is defined as jij,{x) = ql3^3m[{hi — 
ih 2 )D^{hi -\- ih 2 )]- These are precisely equations ( |70| )- 
( [T^ of Section |V] Again, as in the global case, the stan¬ 
dard procedure would be to discretize the continuum ac¬ 
tion such that the covariant derivatives are substituted 
for the standard lattice ones defined in terms of links 
Ui = . However, this would introduce an unnec- 


Pmax — l^min ; L — . (A13) 

2 Pmin 

Note that in this appendix, k refers to physical mo¬ 
mentum and p to lattice momentum. As discussed in 
section m the Higgs EOMs possess a well-known struc¬ 
ture of resonance bands, which can be either of the form 
0 < fc < fc* or of the form fcmin < k < k^. We expect 
these momenta to be physically excited, at least at the 
first stages of the Higgs decay. Therefore, we must have 
a good coverage of this range of momenta. Let us define 
the coefficient ac 


essary complication for describing the term in the 
EOM. Therefore, we proceed again by simply discretizing 
directly the EOM with the correct lattice operators on 
the left-hand side of the equations coming from the dis¬ 
cretization of the lattice gauge invariant action, whereas 
we maintain again the term with the scale factor 
given by a continuous function evaluated at the discrete 
times. The lattice EOMs then look as follows: 

Ag £)-£>+^ 


where we have defined Jn at the lattice point h as 

= ^3m[h*[/oh+o] , (All) 

«o 



Pmin 


The larger is, the better the infrared coverage of the 
resonance band, but the worse the ultraviolet scales are 
captured. In order to probe well the posterior displace¬ 
ment of the spectra to higher momenta when the sys¬ 
tem becomes nonlinear, we need to choose ac judiciously. 
With this idea in mind, we have determined for each g, 
the ac parameter that ensures a good infrared coverage 
without spoiling the ultraviolet part of the spectra. For 
simulations with N = 128 points, we have fixed ac typi- 
’cally within the range 4 < ac ^ H- 

We show in Fig. two particular spectra obtained 
(A10)from simulations of the global scenarios at two differ¬ 
ent times, for different (A, ac) parameters. Apart from 
a better or worse coverage of the ultraviolet or infrared 
regimes, the main physical results are well captured in all 
simulations, and are also consistent between them. The 
same consistency is observed in the Abelian-Higgs simu¬ 
lations, making our results robust versus lattice artifacts. 


Ao A+V,-^(A-A+V,-A+A-V,) = ^3m[h*U,h+i] 

3 

^A-A+V) = Jn , 


and the lattice covariant derivatives as D^(j) = 

and D-(j) = 


One needs to check that for all times, the discrete 
Gauss law (All) is conserved. In particular, we require 
that for all times 


Ag 


1 I Aj Aq Ej - Jn\ 


(A12) 


We have checked that this is indeed the case. In partic¬ 
ular, we find that depending on the simulation, at the 


Appendix B: Initial conditions 

In this section, we discuss in more detail the initial 
conditions of our lattice simulations for both the global 
and the Abelian-Higgs models. As has already been men¬ 
tioned in the main text, we start our simulations (in both 
the global and Abelian-Higgs models) just after inflation 
ends, which we take as the time z = 0. From z = 0 
to z = Zosc, we keep the gauge bosons deactivated, solv¬ 
ing only the Higgs equation for the homogeneous mode. 
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FIG. 17: Top: The spectra of one scalar field (k = 

k/H,) in the global modeling for two different times, 2 = 173 
and 2 = 996, and for different sets of {N, Oc). We have taken 
P = 0.01, q — 14, and a scale factor evolving as a RD universe. 
Bottom: The electric spectra k^\£k,\^ in the Abelian-Higgs 
modeling for the time 2 = 615, for different (X, q^), and for 
P = 0.01, q = 101 and a RD universe. 


Therefore, at 2 = 2 osc, we have h(zosc) given by Eq. (16 1 , 
and the rest of scalar/gauge fields set to 0. It is at this 
time that our lattice simulations truly begin, because at 
this moment we put quantum fluctuations over the homo¬ 
geneous modes of both the Higgs and the decay product 
fields. We now explain how these fluctuations are set in 
both models. 


1. Global model 

For sake of clarity, let us come back temporarily to 
physical variables. Let us use /(x) to denote the quan¬ 
tum fluctuations of a field in position space, and fk its 
Fourier transform, defined as 

At 2 = 2 osc we set, over the homogeneous mode of the 
different fields, a spectrum of quantum fluctuations cor¬ 
responding to the probability distribution of the ground 


state of a scalar field in a FLRW universe 


PilfkMfkl 


2|M - 
(l/fcP) 


ifkp 


where we have 





(B2) 


(B3) 


Here, tOk^osc = \/k'^ + a^c^osc frequency of the 

field at the time 2 osc, and TOosc is the mass at this same 
time, Wosj. = {d‘^V/df^){zosc) with V the potential. The 
mode fk also contains an arbitrary random constant 
phase VArg(/fc). To maintain isotropy properties, we add 
both left-moving and right-moving waves, so that we take 


fk = fk.i + fk,r = , (B4) 

fk = iujka{fk,i - fk,r) - Ufk , (B5) 


where 61 and 62 are constants with 9i G [0, 27r). In the 
discrete lattice, we set the fluctuations in momentum 
space so that, from lat tice point to lattice point, \fk\ 
varies according to Eq. (B2|, and the phases 6*1 and 62 
vary randomly within the interval 6 i G [0, 27r). 

From the properties of the Lame equation discussed in 
Section m we know that depending on the value of the 
resonance parameter q = ( 7 ^/( 4 A), the Higgs equation 
has a certain structure of resonance bands. As we see 
in Fig. the most infrared band is always the one with 
the greatest Floquet index, and we hence expect that the 
Higgs decay will be dominated by this band, at least at 
initial times. It has a maximum at a given momentum, 
which we call fc^ax- Therefore, this allows us to set a 
cutoff to the probability spectrum (B2|, such that for 
k > \fk\ = 0. As it should be, we have confirmed 

that changing this cutoff within a wide range of values 
does not significantly modify our results. 


2. Abelian-Higgs model: Gauss conservation law 


We now discuss how we set the initial quantum fluc¬ 
tuations in the Abelian-Higgs model. Caution must be 
taken in this step, because as we wi ll see, we must ensure 
that the Gauss condition [Eq. (A12) in the discrete] holds 
at the beginning of the simulations. 

Let us come back to natural variables. In this sec¬ 
tion, we define hj{z, 2 osc) with j = 1 , 2 to be the fluctua¬ 
tions of the two components of the conformally rescaled 
Higgs held at the time 2 osc- Let us also define hj{k) = 
hj(k,Zosc) to be their corresponding Fourier transforms. 
Following very closely our discussion of the initial con¬ 
ditions in the global modeling, we impose, over the two 
components of the Higgs, the spectra 


hfk) = 

h2{k) = 

V2 


(B 6 ) 
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Here, \hi\ and |/i 2 | are quantities that change, from 
point to point of the lattice in momentum space, accord¬ 
ing to the probability distribution function 

217? -I - 

= ^d\h,\ , (B7) 


(j = 1 , 2 ) where 




A 


2i?|/32a;, ' 


(B 8 ) 


The frequency is ujj = 
natural momentum, and the natural masses defined as 


rrij, with k = k/H^ the 


— (3^1osc + ^2 osc)/ 5 “ 3^1osc/3 J 

^2 = (^losc + 37l2osc)/3 = ^losci^ ) 


(B9) 


where Tiiosc = hi{zosc) and 7 i 2 osc = h 2 {zosc)- For the 
last equality, we have used that, for the initial conditions 
of the Higgs homogeneous mode given in Eq. (771, we 


have 7 i 2 osc = 0. From (B 6 ), the fluctuations of the Higgs 
derivatives are 

Kik) = - e*®") , 


Kik) = ^*^^2 




(BIO) 


Also, the four different phases vary, in momentum 
space, from lattice point to lattice point. These phases 
would vary in principle randomly within the interval 
9i e [0,27r), but as we are working also with gauge 
bosons, we need to preserve the Gauss law initially. Due 
to this, we thus may need to impose one simple constraint 
to the phases. 

Let us discuss this in more detail. As mentioned before, 
we must ensure at the initial time the Gauss law 


dzV' = joiz) 


(Bll) 


with jo{z) = q/3'^'3m[{hi —ih 2 ){h[ -|-i 7 i 2 )]. Therefore, the 
quantum fluctuations we impose on the gauge fields at 
^osc must preserve this condition. Let us write the Gauss 


law (Bll) in momentum space as 


V-{k,Zosc) = ’i-j^joik) , 

F/(0,Zosc) = 0, 


(B12) 


where jo{k) is the Fourier transform of jo{z) at the time 
2 : = Zosc, and Pmin is the minimum momentum of the 
lattice. This allows us to set fluctuations to the gauge 
fields in the following way: First, for a given lattice point 
in momentum space, we produce the Higgs fluctuations 
according to Eqs. (B 6 ) and (BIO). With these Higgs 


fluctuations, we obtain the correspondent fluctuations of 
jo{z) and its corresponding Fourier trans form joik). Fi¬ 
nally, we fix V-{k,Zosc) according to Eq. (B12). We have 


then obtained a spectrum of initial gauge fluctuations. 

However, in order for this procedure to be valid, we 
must ensure that our current jo{k) does not possess a 
zero mode, i.e. jo(k = 0) = 0. This requirement can 
be clearly seen in Eq. (B12). This is equivalent to say¬ 


ing that there must not be a total electric charge in our 
lattice box. However, from the spectrum of Higgs fluctu¬ 
ations described above, we obtain 


jo{k = 0 ) = / d^zjo{z) = 


(B13) 


= d^kmc[hi{k)h'2{k) - h[{k)h2{k)] 


with 


D\z[hi{k)h2{k) - h2{k)h[{k)] = \hi\\h2\qP‘^ x 


cos I 


^3 + ^4~^l~^2 


) X [a ;2 sin (- 
—wi sin ) cos ( 


2 

O 4 — O 3 
2 


) cos 
)] ■ 


(^) 


(B14) 


This quantity is not zero in general. There does not seem 
to be a particular reason why we should have a total elec¬ 
tric charge in our box, so we should find a way of mak¬ 
ing Eq. (B14) null. We have found two different ways of 


modifying slightly the initial quantum fluctuations of the 
Higgs field to make the integrand of Eq. (B14) zero, which 


do not modify significantly the amplitude of the fluctu¬ 
ations with respect to the approach used in the global 
model. The first one is to impose, at each lattice point, 
the following constraint to the four arbitrary phases of 
the Higgs fluctuations 


04 = 9l + 62 — 03 + TT , 


(B15) 


so that the phases 9i, 02 and 6*3 are randomly generated 
within the interval 0i S [0, 27r), and 6*4 is fixed through 
Eq. (B15|. The second one is to leave the four phases 
totally random, but to perform, at each lattice point, the 
following shift to the Higgs fluctuations: 


Ti'i ^ h[ + 


Jnh 


011-2 


7i? -I- Tio 




hi 


(B16) 


where Jq = (1/A7^) is a sum over all lat¬ 

tice points. This shift eliminates by hand the zero mode 
of the current. One can easily confirm that both meth¬ 
ods make zero the integrand of Eq. (B14). In practice, 
we have confirmed that both methods produce almost 
identical results. This is normal, as in order to trust our 
lattice simulations, the way in which we set the initial 
fluctuations must not play any relevant role, as long as 
their amplitude does not significantly change. 
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